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ABSTRACT 



We introduce an improved code for simulations of star clusters, called MOCCA. 
It combines the Monte Carlo method for star cluster evolution and the Fewbody code 
to perform scattering experiments. The Fewbody was added in order to track more 
precisely dynamical interactions between objects which can lead to creations of various 
exotic objects observed in the star clusters, like Blue Stragglers Stars (BSS). The 
MOCCA code is currently one of the most advanced codes for simulating real size 
star clusters. It follows the star cluster evolution closely to N-body codes but is much 
faster. We show that the MOCCA code is able to follow the evolution of BSS with 
details. It is a suitable tool to perform full scale evolution of real star clusters and 
detail comparison with observations of exotic star cluster objects like BSS. 

This paper is the first one of the series of papers about properties of BSS in star 
clusters. This type of stars is particularly interesting today, because by studying them 
one can get important constrains on a link between the stellar and dynamical evolution 
of star clusters. We discuss here first results concerning BSS for an arbitrary chosen 
test model. We investigate properties of BSS which characterize different channels of 
formation like masses, semi-major axes, eccentricities, and orbital periods. We show 
how BSS from different channels change their types, and discuss initial and final 
positions of BSS, their bimodal distribution in the star cluster, lifetimes and more. 

Key words: 
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1 INTRODUCTION about statistics of BSS. They are particularly interesting 

today, because by studying these types of objects, one can 
This is the first paper of the series of papers about Simula- , . . .. , , . , „ , , 

get important constrains on a link between stellar and dy- 
tions of star clusters, usme improved version of the Monte . , , ,. r , , , 

f namical evolution of the star clusters. Star clusters are very 

Carlo code, called MOCCA, which stands for MOnte Carlo „ . , . c ,. , 

' efficient environments for creating exotic objects like BSS. 

Cluster simulAtor. The MOCCA code is currently one of , . it , , . , , . . r 

t l i- By studying them one can reveal the dynamical history of a 



the most advanced codes which is able to simulate real size 
star clusters and at the same time, it allows to have a full 
dynamical history of the evolution of all stars in the system. 



cluster and the role of dynamics on the stellar evolution. BSS 

statistics can also provide some constrains for initial binary 

properties. Thus, numerical codes which are able to simulate 
It follows the star cluster evolution closely to N-body codes , , , , . r . , , , 

J the complete evolution of star clusters and at the same time 

but is much faster. We show that the MOCCA code is able P „ , . , . , , , . . „ 

follow movements, dynamical interactions and evolution of 
to follow the evolution of exotic objects, m this case Blue , „ , . .„ , „ ,, r , . 

. any stellar objects, are very significant. Results of such sim- 



ulations, by comparing with observational data, can verify 



Stragglers Stars (BSS), with details and thus it is a suitable 

tool to perform meaningful analysis and comparisons with 

• r i i • i „. : n many theories. 

observ ations of exotic star cluster objects. Chat teriee et al.l 

|2010D are among the first who showed that the Monte Carlo BSS are defined as stars which are brighter and bluer 

method can be a very good tool for studying BSS. (hotter) than the main sequence turn-off point. These stars 

The main subject of this paper concerns first results he along an extension of the main sequence (MS) in the 

Color-Magnitude Diagram (CMD) and appear to be reju- 
venated stellar population. BSS are on the place in CMD 
* E-mail: ahypki@camk.edu.pl where they should already evolve away from the main se- 
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quence. Their mass is larger that the turn-off mass, which 
suggests some stellar merger or a mass tran sfer scenario for 
their creation. They were first discovered bv lSandagd (|l953l ) 
in M3 and later observations s howed that BSS ar e present 
essentially in all star clusters. (|Piotto et al.l I2004T ) counted 
3000 BSS in 56 different size clusters. 

Currently, there are two main scenarios considered as 
possible formation mechanisms for BSS. The first scenario 
is a mass transfer between binary compani ons which can 
lead to the coalescen ce of the binary system (|McCrealll96i 
IZinn fc Searie1ll976r h The second leading scen ario for cre- 
ating BSS is a physical collision between stars l|Hills fc Davl 
1 19761 ) . However, the exact nature of channe ls of formation of 
these objects is still unclear. According to iFusi Pecci et all 
<|l992h . different environments could be responsible for dif- 
ferent origins of BSS. In globular clusters (GCs) which are 
not dense, BSS could form as evolutionary mergers of pri- 
mordial binaries, and in high density GCs, BSS could form 
from dynamical interactions, particularly from interactions 
involving binaries. 

Relative efficiency of these two main formation channels 
is still unknown. Though, it is believed, that they act with 
different efficiency according to the cluster structural pa- 
rameters (|Fusi Pecci et al.lll992T l and additionally they can 
work simultaneously in different radial pa rts of a star cluster 
ijFerraro et al.lll997i : iMapelli et all 120061 ). Particularly, the 
number of BSS formed in th e cluster does not correlate with 
the pr edicted collision rate l|Piotto et al.l |2004| ; iLeigh et al.l 

200830)- 

This is one of the reasons why it is believed that 
mass transfer mechanism is more imp ortant in creation o f 
BSS instead of collision between stars l|Knigge et al.l l2009). 
Unfortunately, there is still no simple observational distinc- 
tion between BSS formation through mass transfer or colli- 
sion between stars. One of the first attemp t s to c larify this 
issue is the approach of lFerraro fc Lanzonil (|2009l ) . who ob- 
served a significant depletion of C and O suggesting mass 
transfer mechanism for creating some BSS sub-population 
in 47 Tuc. According to iDavies et al.l (|2004T ) primordial bi- 
naries with BSS are vulnerable to exchange encounters in 
the crowded environments of star clusters. Low-mass compo- 
nents are replaced by more massive single stars. The authors 
claim that these encounters tend to reduce the number of bi- 
naries containing primaries with masses close to the present 
turn-off mass. Thus, the population of primordial BSS is 
reduced in more massive star clusters. 

iFerraro eT al. (2003) defined the BSS specific frequency 
as the number of BSS, normalized to the number of the 
horizontal branch stars. They examined 6 GCs and found 
that BSS specific frequency varies from 0.07 to 0.92, and 
it does not depend on central density, total mass and ve- 
locity dispersion. What is surprising, they found the largest 
BSS specific frequencies for clusters with the lowest central 
density (NGC 288) and the highest central density (M80). 
IFerraro et all (|2003T ) claim that these two kinds of BSS for- 
mation processes, mass transfer and mergers, can have com- 
parable efficiency i n producing BSS in their respective typ- 
ical environments. ISollima et al.l (|2008T ) found a strong cor- 
relation between BSS specific frequency and linear combi- 
nation of binary fraction (^Mn) and velocity dispersion (a v ) 
€,bin + cea v where a = —4.62. This indicates that, for a given 
binary fraction, BSS specific frequency decreases with in- 
creasing velocity dispersion. Small cluster velocity dispersion 



corresponds to a lower binding energy limit between soft and 
hard binaries (to a larger fraction of hard binaries). Since 
the natural evolution of hard bina ries leads to increase of 
their binding energy (Heggic 1975), low velocity dispersion 
GCs should host a larger fraction of hard binaries, which are 
able to both survive possible stellar encounters, and activate 
mass-tra nsfer and /or mergin g processes between the com- 
panions (|Sollima et al.l 120081 ). Therefore, more BSS formed 
by the evolution of primordial binaries are expected to form 
in lower velocity dispersion GCs. 

ISoilima"eT al. (2008) tested 13 low-density GCs for cor- 
relations between specific frequency of BSS and cluster pa- 
rameters like binary fraction, total magnitude, age, central 
velocity dispersion, metallicity, cluster central density, half- 
mass relaxation time, half-mass radius, stellar collision rate, 
concentration, and cluster evaporation rate. BSS specific fre- 
quency was defined as ratio between estimated BSS num- 
ber and MS number. MS were chosen, instead of horizontal 
branch (HB) or red giant branch (RGB) stars, because of 
their abundance in all clusters and its completeness. They 
found the strongest correlation between number of BSS and 
binary fraction. It suggests that the primordial binaries frac- 
tion is one of the most important factor for producing BSS. 
Additionally, noticeable correlation exists with the absolute 
magnitude and anticorrelation with the cluster age and cen- 
tral velocity dispersion. The age estimates are uncertain and 
span a narrow range, so one has to be careful while mak- 
ing some generalizations. However, if such anticorrelation 
with the cluster ages would be confirmed in the future, it 
could suggest that binaries disruptions in cores of GCs be- 
come more efficient with time, which would in consequence 
reduce the fraction of binaries and also BSS in the core. 
ISollima et al.l (|2008T ) suggest that strong correlation between 
number of BSS and binary fraction is a result of formation 
channel of BSS as the unperturbed evolution of primordial 
binary systems. They found no correlations for central den- 
sity, concentration, stellar collision rate and half-mass relax- 
ation time. This indicates that the collisional channel of BSS 
formation has very small efficiency in low-density GCs. 

Radial distribution of BSS in many clusters is bimodal. 
First discoveries of bimodal dist ributio ns were done f or the 
M3 bv lFerraro et al. | l|l993l , ll997ft and bv lZaggia etail (|l997t ) 
for M55. BSS radial distribution for M3 cluster clearly shows 
the maximum at the center of the cluster, clear-cut dip in 
the intermediate region and again rise of BSS in the outer 
region of the cluster (but lower than the central value). Bi- 
modal distribution of BSS was later shown by other au - 
thors for ot her clusters like 47 Tuc (Fcrraro et al. 2004), 
NG C 6752 JSabbi et al l l2004h, M55 jLanzoni et al l I2007T ) , 



M5 dWarren et all 20061 ; lLanzoni et al.ll2007i ), and others. 
ISigurdsson et al. n|l994h suggested that BSS were formed by 
direct collisions in the center of the cluster and then ejected 
to the outer part of the system as a result of a dynamical 
interaction. Ejected BSS would afterwards moved back to 
the center of the cluster because of the mass segregation, 
which leads to the increase of the number of BSS in the cen- 
ter of the system. If the dynamical interaction is energetic 
enough then BSS would stay outside of the cluster for longer 
time and this could be the reason why there is a higher rate 
of BSS in the outer part of the cluster - second peak of 
BSS in bimodal distribution. Later, bimodal distribution of 
BSS in the cluster was explained differently bv lFerraro et al.l 
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| 19971 ). They showed it is a result of different processes form- 
ing BSS in different parts of the cluster - mass transfer for 
the outer BSS and stellar collisions leading to mergers for 
BSS in the c enter of the cluster. F urthe rmore jMapelli et all 
|2004 I2006T I and lLanzoni et all |2007fi by performing nu- 
merical simulations showed that bimodal distribution in the 
cluster cannot be explained only by collisional scenario in 
which BSS are created in the center of the cluster and some 
of them are ejected to the outer part of the system. This pro- 
cess is believed to be not efficient enough, and ~ 20 — 40% 
of BSS have to be created in the peripherals in order to get 
the required number of BSS for the cluster, ft is believed, 
that in the outer part of a star cluster, binaries can start 
mass transfer in isolation without suffering from energetic 
dynamical interactions with field stars. Even if one can ob- 
serve bimodal distribution of BSS for many clusters, one 
can not generalize this feature. There are known clusters for 
which radial distributions are even flat, like for NGC 2419 
l|Dalessandro et alj|200i ; IContreras Ramos et ai1l2012l) . 

iFerraro et all (|2006h gave first results of currently being 
performed extensive surveys of chemical composition of BSS 
for some selected GCs. They examined 43 BSS in 47 Tuc 
and found the first evidence that some sub-population of 
these BSS have significant depletion of C and O with re- 
spect to the normal cluster stars. They argue that this is 
caused by CNO burning products on BSS surface, coming 
from the core of a deeply peeled primary star. This scenario 
is expected for the case of mass transfer formation mecha- 
nism and coul d be the first direct proof of this formation 
process. Later, iFossati et all (|201Ch attempted to develop a 
formation scenario for HD 73666, a known BSS from the 
Praesepe cluster, and showed that abundance of CNO is 
consistent with a collisional formation. However, they were 
unable to determine whether HD 73666 is a product of a 
collision between two stars, components of a binary or be- 
tween binary systems. Further studies of these phenomena 
could create some statistics on how efficient this mechanism 
coul d be in produ c ing B S S . 

IFerraro et ail (|2009l ) reported two distinct sequences of 
BSS in the globular cluster M30. These two groups are 
clearly separate in the C MP and nearly parallel to each 
other (jFerraro et al.ll2009l . Fig. 1). The first BSS sequence 
was accurately rep roduced by the collisional isochrones 
|Ferraro et alj 120091 Fig. 4, blue points). The second BSS 
sequence well corresponds to the ZAMS shifted by 0.75 
mag, marking the position of the low-luminosity boundary 
predicted for a pop ulation of mass-transfer binary systems 
|Ferraro et ai1l2009l, Fig . 4, red points). 

iKnigge et al.l (|2009h focused on BSS in cores of star 
clusters, because in these regions collisions between stars 
should be the most frequent. They used existing data from 
a large set of HST-based CMDs and confirmed that there is 
no global correlation between the observed core BSS number 
and the collision rate (different core densities have different 
predicted collision rates and it does not correlate with the 
number of BSS). However, there is a significant correlation 
if one wo uld restrict thi s relat ion to the clusters with dense 
cores (see IKnigge et all |2009h black points in Fig. 1). The 
second relation which was tested by this group concerns the 
binary fraction in the core. If most of BSS are formed in 
binaries, the number of BSS should scale with the binary 
fraction simply as Nbss oc fbinM core , where fbin is the 



binary fraction in the core, and M core is the total stellar 
mass contained in the core. Indeed, they found clear corre- 
lation between the number of BSS and core masses of the 
clusters, as it is expected for the scenario where most BSS 
originate from binaries (see IKnigge et alj (|2009l . Fig. 2)). 
They interpret this result as a strong evidence that more 
BSS originates from binaries instead of collisions between 
stars. They found that dependence Nbss oc M c s ore can be 
estimated with 8 ~ 0.4 — 0.5. Furthermore, they estimated 
power law co r relatio n fa„ oc M~ a ®f 5 based on the data from 
iMilone et al.l ((2008) who described global parameters for 35 
clusters spanning a wide range of density and other dynam- 
ical star cluster parameters. Those two estimates combined 
together shows that the number of BSS found in the cores of 
globular clusters scales roughly as Nbss oc fbinMcore, just 
as expected if most core BSS are formed in binary systems 
|Knigge et alj|2009h . 

BSS are being found in the halo and in the bulge of 
the Galaxy (iBragaglia et al. 200? ; Fuhrmann et al.l 1201 ll ; 
IClarkson et al.ll201ll ). iTillich et alJ(|201Ch found that a star 
SDSSJ130005. 62+042201.6 (J1300+0422 for short) is a BSS 
from the halo and has radial velocity of about 504.6 ± 
5 km/s. With Galactic rest-frame velocity of about 467 
km/s, J1300+0422 travels faster than any known blue strag- 
gler, but still is bound to the Galaxy. 

Recently, iGeller fc Mathieul |201ll ) reported that BSS 
in long period binaries in an old (7 Gyr) open cluster, 
NGC 188, have companions with masses of about half of the 
solar mass, which is a surprisingly narrow mass distribution. 
This rules out a collisional origin for these long period BSS, 
because otherwise, for collision hypothesis there would be 
significantly more companions with higher masses. The data 
is consistent with a mass transfer origin for the long-period 
blue straggler binaries in NGC 188, in which the compan- 
ions would be white dw arfs of about half of a solar mass 
|Geller fc Mathieull201lh . 

At the end of this section we would like to note that 
in the literature, terms collision and merger are used dif- 
ferently by different authors. In this work the term collision 
is defined as a physical collision between at least two stars 
during a dynamical interaction, while the term merger is 
defined as a coalescence between stars from one binary as a 
result of stellar evolution. 

This paper is organized as follows. In the Sec.[2]there is 
described the MOCCA code, its advantages and drawbacks 
in comparison to N-body codes, its design and interface be- 
tween different internal parts of the code. In the Sec.[3]there 
is a description of the initial conditions for the test simula- 
tion, used in this paper, which shows the ability of the code 
to follow the evolution of exotic objects (BSS in this case). 
Sec. [4] contains first results of the simulation and physical 
interpretation of properties of BSS for different channels of 
formation. We will discuss BSS masses, orbital periods, ec- 
centricities, BSS type changes and possible induced mass 
transfer which could lead to BSS formation. Moreover, we 
will discuss BSS global properties, like their initial and fi- 
nal positions in the star cluster, BSS bimodal distribution 
and their lifetimes. Final Sec. [S] summarizes our findings, 
presents discussion about channels of formation of BSS, and 
highlights the future plans for the MOCCA code. 
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2 DESCRIPTION OF THE CODE 

We give here a description of internal parts of the MOCCA 
code and the interface which connects them into one pack- 
age. In the first subsection there is described the old ver- 
sion of the code, which is basically Monte Carlo method for 
star clusters simulations together with dynamical interac- 
tions based on cross-sections, and cod e for stellar evolution . 
Then, we describe the Fewbody code (|Freeeau et al.ll2004h . 
which was added to deal with the dynamical interactions 
between binaries and single stars or between two binaries. 
Next, we give the technical reference of the interface which 
connects these two codes. In the last subsection, there is a 
brief summary about the performance of the code. 

This paper introduces for the first time the MOCCA 
code. Thus, in this section there is a detailed technical 
description of its internal parts, particularly the interface 
which connects them. This description is important to un- 
derstand the design of the code and to be able to develop 
new procedures into the MOCCA code. 



2.1 Old version of the code 

The old version of the code refers to the code which uses 
Monte Carlo approach to describe the relaxation processes 
for star clusters but with some additional features which 
have been developed over the years by Giersz and h is col- 
laborators |Giersj|l998l . l200ll , l2006l ; lGiersz et afll2008h . The 
code for Monte Carlo method itself is in turn based on the 
orbit-averag ed Mon t e Car lo method developed in the early 
se venties by Henon ( 197lh and then substantially improved 
bv lStodolkiewic J ( 19861 ). Additionally, the old version of the 
code is equipped with procedures for performing three- and 
four-body dynamical interactions, based on a cross-section 
approach, involving primordial binaries and binaries formed 
dynamically. It takes into account the Galactic tide dur- 
ing escape process es, according to the theory proposed by 
iBaumgardtl l|200ll ). and it has implemented procedures to 
perform th e internal stellar evolutio n of both single and bi- 
nary stars (|Hurlev et al.ll2000l . 120021). Details about th e old 
version of the code one can find in (|Giersz et al.ll201ll . and 
references within). 

The old version of the code had several shortcom- 
ings. Dynamical interactions between objects were calcu- 
lated analytically, using cross-sections. The problem is that 
cross-sections are not well known in the case of unequal 
masses and complicated resonant interactions. Additionally, 
for cross-sections possibility of stellar collisions are excluded. 
This is one of the most important shortcomings of the old 
version of the Monte Carlo code. The Fewbody code fixes this 
problem and in the current version of the MOCCA code 
there are possible all kinds of outcomes for dynamical in- 
teractions. Moreover, the Fewbody allows to have binaries 
hardening and softening with respect to only hardening in 
the old version of the code. This is a very important im- 
provement which opens a lot of new interaction channels 
and allows to precisely track changes of energy of binaries. 
Details about the advantages of using the Fewbody one can 
find in the Sec. H2] 

The second main shortcoming of the old version of the 
Monte Carlo code concerned escaping stars from a tidal lim- 
ited star cluster. The escape rate was scaling with N 3 ^ 4 to 



deal with backscattering process (|Baumgardtll200ll ) but a 
star or a binary was removed from a star cluster imme- 
diately. To describe more realistically the escape process 
we im plemented procedures based on iFukushige fc Heggid 
(2000). The escape is not anymore immediate but stars need 
some time to escape from a star cluster. It can significantly 
delay the time of the escape of a star from the star clus- 
ter. Details one ca n find in the next paper of this series 
(|Giersz et alj|201ll ). 



2.2 The Fewbody code 

The Fewbody (jFregeau et al.l I2004T ) is a software package 
for performing small-N scattering experiments. It uses the 
8-th order Runge-Kutta Prince-Dormand integrator to ad- 
vance the particles' positions. There is a possibility to en- 
able the_Jull_p_airwise_JK-S regularization in the simulation 
too (|Aarseth fc Zardll974l ). The Fewbody code detects stable 
hierarchical systems and isolates unperturbed hierarchies to 
increase dramatically overall performance. Hierarchies and 
internal data structures of stars are stored in binary trees 
which means that each bound object can have only two child 
objects (t he simplest hierarchy is a b inary star). The Few- 
body uses iMardling fc Aarsethl (|200lh stability criterion to 
assess the stability of hierarchies at each level and interrupts 
calculations if all bound objects are considered as stable. 
This code can handle dynamical interactions between any 
arbitrary number of stars and understands arbitrary com- 
plicated hi erarchies. Full details about the Fewbody code one 
can find in lFregeau et al.l (2004). 

The Fewbody, in the present version of the MOCCA 
code, is used to perform binary-single and binary-binary 
interactions. Incorporating it for performing scattering dy- 
namical interactions, brings additional desired features. The 
MOCCA code runs independently small N-body scattering 
experiments like binary-binary interactions, using the Few- 
body, and afterwards results are written back to the MOCCA 
data structures. 

The Fewbody in the MOCCA code allows to follow the 
dynamical evolution of binary objects and higher hierarchies 
just like in N-body simulations, but still giving opportunity 
to run simulations several orders of magnitudes faster than 
the N-body codes (see Sec. 12. 5p . The Fewbody allows the 
MOCCA code to perform interactions between stars and bi- 
naries more realistically than using analytical formulas like 
in the old version. In the beginning of the clusters simu- 
lations the most important process is the stellar evolution 
but after some time, in many clusters, dynamical interac- 
tions between stars start to play a huge role in the overall 
clusters evolution. Thus, dynamical interactions can be very 
significant not only for the global cluster evolution but also 
if one has to study creation and evolution of many different, 
exotic objects like compact binaries, black holes, BSS and 
more. Thus, the Fewbody is so much needed to have more 
realistic simulations. 

The channels of formation of BSS are the main subject 
investigated in this paper. The Fewbody allows the MOCCA 
code to have channels of formation essentially the same as 
in N-body simulations. It is possible to have any outcome 
from interactions between stars and binaries, like single or 
multiple mergers, exchanges and disruptions, both in simple 
dynamical interactions and for resonant ones. Resonant in- 
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teractions occur when an incoming star is temporarily bound 
to the binary during the scattering experiment. Eventually 
one star has to escape, but before it happens, some complex 
and chaotic stars movements can occur. In the old version 
of the code it was not possible to get all complex outcomes 
from dynamical interactions. In comparison to the old ver- 
sion of the code, we expect to have similar global properties 
of star clusters, but parameters which strongly depend on 
dynamical interactions (like number of BSS or spatial dis- 
tribution of BSS) should be closer to observations or results 
coming from pure N-body simulations. 

Input parameters for the Fewbody are stars' and bina- 
ries' parameters such as masses rrii, radii Ri, semi-major 
axes ai, eccentricities ej, and some global values which char- 
acterize a dynamical interaction like impact parameter b, 
relative velocity at infinity Voo, and technical parameters 
like tidal perturbation, maximum time for computations in 
seconds tcpu, dynamical time in years td yn , or KS regular- 
ization - a coordinate transformation that removes all singu- 
larities from the N-body equations, making the integration 
of close approaches much more accurate. 

Impact parameter, b, is divided by a sum of the semi- 
major axes of all interacting binaries. Relative velocity, v r , 
is the relative velocity at the infinity between interacting 
distinct bound objects in terms of the critical velocity, v c , 
defined in a way that the total energy of the binary-single 
or binary-binary system is zero. If Voo > v c the total energy 
of the system is positive and it is possible that each object 
will leave the system unbound from any other with positive 
velocity at infinity. If Voo < v c the total energy is negative 
and the encounters are likely to be resonant, with all stars 
involved remaining in a small volume for many dynamical 
times. 

Tidal perturbation determines when analytic formulas 
and direct integration procedures are used. Each numeri- 
cal integration is started at the point at which the tidal 
perturbation on a binary in the system reaches some spec- 
ified value 5 (S = 10" 5 is the default value in the MOCCA 
code). Tidal perturbation is defined as Fud/F re i, where Fud 
is the tidal force at the apocentre, and F re i is the relati ve 
force at the apocentre (for details see lFregeau et af] l|2004h ). 
The same mechanism is used internally by the Fewbody to 
speed up integration between stars and bound objects. The 
force is not calculated between all stars but between objects 
which do not break this tidal perturbation threshold. It is a 
quite important parameter because smaller values of S yield 
better energy conservation but increase the computational 
time - more integration steps are calculated with numerical 
integrator rather than with analytical equat ions (for more 
inform ation about the hierarchy isolation see lFregeau et all 
<l2004h ). 

Maximum time of computations in seconds for each 
Fewbody scattering experiment is set by parameter tcpu- 
After this time the interaction is forced to stop. It is pos- 
sible that stars are still close to each other (tidal pertur- 
bation is still > 8). Thus, this parameter has to be cho- 
sen carefully. It can not be too short because many more 
dynamical interactions would not be completed according 
to stopping conditions (described later in this section). By 
experimenting and calculating how many interactions were 
not completed, tcpu = 10 s was chosen as an optimal value. 
For maximum time 10 seconds for a simulation with 100k 



objects there were only few dozens of not completed en- 
counters (from total 4057 interactions) - it seems to be a 
reasonable value. It is also worth to mention, that it is hard 
to predict in advance how many interaction steps each scat- 
tering experiment would take, so a better solution is to use 
tcpu- Interrupted interactions resulted in creating unsta- 
ble or stable triples or quadruples. These objects were ar- 
tificially disrupted to binaries and single stars because the 
Monte Carlo part of the MOCCA code is currently unable to 
handle complex hierarchies. However, even if those objects 
were manually disrupted, the binding energy of triples and 
quadruples were insignificant in comparison to the average 
binding energy of binaries in the system. Thus manual dis- 
ruption most probably had no significant influence on the 
overall cluster simulation. 

KS regularization for our simulations is disabled by de- 
fault. Regularization transforms coordinates of stars remov- 
ing all singularities from N-body equations. It allows to inte- 
grate close approaches and even collision orbits much more 
accurate. But using the KS regularization requires addi- 
tional efforts to detect physical collisions, because pericenter 
is not necessarily resolved by the integrator. With enabled 
KS regularization there were only a few mergers per sim- 
ulation. Thus, we decided to switch it off completely. For 
more informat i on ab out the Fewbody input parameters see 
iFregeau et ail (|2004T l. 

When describing the Fewbody one has to realize what 
are the Fewbody stopping conditions, because it is crucial 
to find optimal parameters for the MOCCA code to have 
a good error conservation with decrease in performance as 
low as possible. Better energy conservation is for smaller 
S, but the Fewbody dynamical interactions are calculated 
significantly longer. The Fewbody uses several criteria to au- 
tomatically terminate the integration of the scattering en- 
counter. In general calculations are interrupted when there 
is no chance to bound objects to interact with each other, 
and bound object will not evolve int ernally. Stopping con - 
ditions are described in details in the Freg eau et alj l| 20041) . 

2.3 Interface design 

Writing an interface between Monte Carlo code and the Few- 
body code was a very complicated task. These codes have dif- 
ferent data types storing data, different programming styles, 
and what is the most important, are written in different pro- 
gramming languages. The Fewbody is written in C language 
and the Monte Carlo code is written in Fortran77. 

The approach to describe tree structures in the code, 
without actually having structures in Fortran77, was to use 
the same idea as in relational models for database manage- 
ment. It is a database model based on the fir st-order pred - 
icate logic, first formulated and proposed by ICoddl (|l969l ). 
Each table describes one data type, in other words, one en- 
tity with some finite set of attributes. Each attribute corre- 
sponds to a column in such a database and each data type 
corresponds to a table. There are many advantages by or- 
ganizing data in such a way. Data types are easy to extend 
in both directions - higher hierarchies (triples, quadruples 
etc.) and more attributes describing objects (like mass, ra- 
dius, luminosity for stars, or semi-major axis, eccentricity for 
binaries). If one has to add some more attributes or more 
objects for an encounter, there is no need to change func- 
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# Data Parent Childl 



Child2 New Merger 
parent id 



-1 
-2 
-3 
-1 
1 
2 
3 








-1 

-3 
-1 







-2 
-4 
-5 



Table 1. Table hierarchy holds information about the structure 
of bound objects 



# a [AU] 



Binary id New binary id 



1.60 
0.98 
1.80 



0.46 
0.31 
0.95 



119390 
825260 




119390 


119390 



Table 2. Table binary holds physical properties of the bound 
objects like semi-major axes (a) and eccentricities (e). The sim- 
plest bound object is a binary. Values of semi-major axes and 
eccentricities are exemplary. 



# 


Mass [M Q ] 


Radius [-Rq] 


Star id 


1 


0.76 


0.70 


119390 


2 


0.75 


0.69 


119391 


3 


0.98 


0.88 


825260 


4 


0.60 


0.53 


825270 


5 


1.74 


4.74 


119391 



Table 3. Table single holds information about physical proper- 
ties of the stars. Masses, radii and star IDs' are exemplary. 



tion declarations, but only extend what is already written 
to deal with new attributes. Especially, there is no need to 
change function declarations by adding or removing input 
parameters. Using such data organization one can have easy 
and transparent system to exchange data between very dif- 
ferent codes. It is also worth to notice that this is just a 
concept of storing data. It does not involve using any un- 
usual Fortran77 language statements. The interface simply 
uses plain arrays of double precision and integer numbers 
to store data, runs the Fewbody, and rewrites results back 
to Monte Carlo structures. From now on, the term 'table' 
means actually just simple 2-dimensional array. 

There are defined only 3 tables to store data: hierarchy, 
binary, and single. This is enough to describe any complex 
hierarchical system with any number of stars in it. 

Detailed description of these tables will be presented 
by the following example. Let us consider an interaction be- 
tween two binary stars. Let us assume that the result of 
such interaction is one binary with one merger inside and 
one single stars (outside the binary). 

The data in the table hierarchy for such encounter is 
shown in Tab. [I] It contains six columns: data, parent, childl, 
child2, new parent and merger id. In this table there are 
stored relations between stars and it contains information 
on which star belongs to which binary. Additionally, each 
row in this table describes what are the children of a specific 



object and which physical properties from two other tables 
[binary or single) describe this object. 

In order to describe a binary there are needed three 
rows in the hierarchy table. Column data contains row ID 
of the data for a single star or for a binary. In our example 
in Tab. [T] for the two first rows we have in the column data 
indices -1 and -2. Negative values mean that these are dis- 
tinct stars and their absolute values point to the first and the 
second row in the single table (Tab. [3}. And if the column 
data contains positive number then it points to the binary 
table (Tab. [21, which contains data of binaries or higher 
hierarchies. Two first rows in the hierarchy table describe 
both single stars of the first binary, there are no children for 
them, and that is why columns childl and child2 contain 0. 

Column parent points to the row in the same table hi- 
erarchy and describes the parent of a given object. Both 
single star and binary can have parents. Of course for the 
simplest case, the binary-single dynamical interaction, only 
single stars (from a binary) can have some indices in the col- 
umn parent. Two first rows in Tab. [1] contain in the column 
parent values 5 which corresponds to the fifth row in the hi- 
erarchy table. It means that these two stars have the same 
parent and they form one binary described by the fifth row 
in this table. This fifth row holds information related to the 
binary itself. Furthermore, in the fifth row of the hierarchy 
table one can see that in the column data there is positive 
index 1. It points to the first row in the binary table (Tab.[2]| 
where all physical properties of such binary are stored, like 
semi-major axis and eccentricity. 

Third and forth row in the hierarchy table describes 
the first and the second star, which are bound together as 
a second binary. This binary containing those two stars is 
described by the sixth row in the hierarchy table. One can 
see that in the sixth row there is the value 2 which is positive 
and thus points to the binary table to the second row. In this 
second row there are stored initial binary properties about 
this second binary taking part in this dynamical interaction. 

The last very important column is merger id. It points 
to the single table in the case of merger. In our example there 
is one merger inside a binary after the interaction. One can 
see that in the column new parent, which describes structure 
after the interaction, in the hierarchy table there are values 
7 for three first rows. It means that those three single stars 
have the same parent. Additionally, in the column merger 
id for the second and the third row, there is the value 5. 
It means that those two stars merged into one star in the 
interaction. Fifth row in the single table describes physical 
properties of this merger. Thus, after the interaction the 
binary consists of two stars where the first star comes from 
the first row {new parent column contains value 7, merger 
id column is equal to 0). The second star is a merger star, 
produced from stars from rows 2 and 3 [merger id column 
is equal to 5). Using these several columns one can describe 
any hierarchical structure needed in the MOCCA code. 

Let us now consider the last two tables (simpler ones). 
The table binary (Tab. [2} is a table (array) which stores in- 
formation about two distinct objects bound gravitationally. 
Such a bound object can be described using attributes like 
semi-major axis and eccentricity. The easiest case is a binary 
star, the more complicated are triples or quadruples. But in 
the case of triple there would be actually two binaries, one 
for the inner binary and the second for the outer binary, 
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so there would be two rows in the binary table needed to 
describe it. 

The table single holds properties of distinct stars 
(Tab. 0J. It contains data like masses, radii, IDs etc. Each 
row in this table corresponds to only one star, merger or 
not, involved in some binary system or not. 

Rows' orders do not matter in any table. In the Tab. fT] 
two first rows are occupied by stars belonging to the first bi- 
nary, the third and forth row describe stars from the second 
binary. However such order is not mandatory. 

The introduced example from tables [1] [2] and [3] shows 
only several columns which are necessary to understand the 
concept behind this interface. In these tables there is stored 
much more information about stars and binaries. For stars 
there are stored additionally velocities from before and after 
interaction. For binaries there are stored additionally bind- 
ing energies, semi-major axes and eccentricities from before 
and after the dynamical interaction. All these pieces of data 
are needed to check how properties of stars and binaries are 
changed after performing one Fewbody dynamical interac- 
tion. All changed properties of the stars and binaries have 
to be saved to the MOCCA structures so that the Monte 
Carlo part of the code will be notified about these changes. 

The source code of the Fewbody is designed very well 
because it is stateless. It does not use any global variables 
(except of course for constant parameters). There was no 
need to worry about global variables when the interface was 
designed, and there was no need to change the Fewbody 
code in order to incorporate it to the MOCCA code. Ev- 
erything what was needed to start the Fewbody was to run 
it with some initial parameters, taken from MOCCA, and 
then read results and propagate them into the MOCCA data 
structures (together with transformation to star cluster co- 
ordinates). Interface design in overall is rather complicated, 
however executing the Fewbody dynamical interaction itself, 
is simple. The Fewbody code design will simplify also the 
parallelization of the code which is planned as one of the 
next new features of the code (discussed later). 

Repeatable results for the simulations with exactly the 
same initial conditions are needed for testing purposes. Pa- 
rameter, which causes that two simulations could calculate 
a little bit differently, is the maximum time for calculations 
in seconds (tcpu)- This value has further implication to the 
whole simulation. On different computers, during those 10 
seconds, integrator in the Fewbody can compute different 
number of integration steps. This will cause, that for ex- 
actly the same input parameters, the same encounter can 
return different results. Usually the difference concerns very 
small changes of such parameters like semi-major axes or 
output velocities. But sometimes the difference can be more 
significant. In the slower CPU during 10 seconds the out- 
come can be for example an unstable triple, but on the faster 
CPU, which is able to calculate more integration steps in 10 
seconds, this unstable triple can be disrupted. Afterwards 
this small difference propagate to the whole cluster. Even a 
tiny difference causes that the whole simulation is not re- 
peatable. It complicates debugging and testing. Statistically 
those changes have no meaning at all, the star cluster global 
parameters are almost the same, but it should be kept in 
mind that two exactly the same simulations, running even 
on the same machine, can give slightly different results. 



2.4 MOCCA code 

The MOCCA code is a package which combines together 
several other codes and allows to perform simulations of 
a real size star clusters. The old version was described in 
Sec. 12.11 Interface which combines the old version of the 
code with Fewbody was described in Sec. 12.31 All these codes 
together create one new package, called the MOCCA code. 

In order to combine Monte Carlo code with Fewbody 
there had to be made some simplifications in the code. As 
was mentioned before, Fewbody can handle any arbitrary 
hierarchy but not Monte Carlo. Therefore, if one of the out- 
comes of the Fewbody is a hierarchical object more compli- 
cated than binary (triple or quadruple), it has to be manu- 
ally disrupted. It has to be done in such a way that overall 
energy has to be conserved. Such cases are only count down 
and written to the log file. However, it is planned to adopt 
the MOCCA code for more complex hierarchies in the fu- 
ture. 

MOCCA code runs each dynamical interaction inde- 
pendently of any other interactions and of the host star 
cluster. Because of that, there are possible dynamical in- 
teractions which create very wide binaries. From Fewbody 
point of view, those objects, even if they are extremely wide, 
are stable and they are not disrupted by the Fewbody code. 
However, in star cluster one has to take into account the 
environment. If semi-major axes of such wide binaries ex- 
ceed many times the mean distance between stars, it is most 
likely that these binaries would not exist at all and would be 
disrupted by encounters right away. Additionally, if binaries 
are very wide then the probability of dynamical interactions 
for them are very high and thus such binaries practically 
always take part in interactions. Furthermore, if semi-major 
axis is very wide, then drawn impact parameter is very large 
and the interaction is just a very distant fly-by (like in re- 
laxation procedures) . Those kinds of dynamical interactions 
are physically unimportant and occupy CPU for no reason. 
Thus, there were implemented procedures to disrupt those 
very wide bi naries in bina ry-single dynamical interactions, 
according to iHeggiel (|l975l . eq. 4.12) probabilities formulas. 



2.5 Code performance 

MOCCA code is much faster than N-body codes for the same 
number of particles in the system. This is simply the con- 
sequence of the Monte Carlo method used in the MOCCA 
code. For the test model, described in details in the Sec. 13. II 
with 120k stars, 20% of primordial binaries, enabled Few- 
body for dynamical interactions, and enabled stellar evolu- 
tion simulation, takes less than 4 hours to complete the cal- 
culations. Simulation with the same initial conditions but 
without using Fewbody for dynamical interactions takes a 
little bit more than 2 hours. The computer used for simula- 
tions had AMD Opteron CPU and 4 GB of memory, so it 
was a simple commodity hardware. N-body codes take weeks 
of computation for similar models wit h even less initial bi- 
nary fractions (e.g. lHurlev et al.l (|2008f )'l. The MOCCA code 
is still under heavy development, it is not optimized yet and 
works only with one CPU. Parallelization is one of the next 
planned features which increase the value of the MOCCA 
code even more. 

The speed of the MOCCA code is its great advantage in 
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Parameter 



Description 



Single stars (N s ) 

Binary stars (iVj,) 

Binary fraction (/(, 

N * ) 
N b +N s I 

Initial model 
IMF of stars 

IMF of binaries 



Total mass (M(0)) 
Observational initial 
core radius (r c i) 
Initial half-light radius 

( r hl) 

Initial tidal radius 
Binary mass ratios 
Binary semi-major 
axes 

Binary eccentricities 



cq. 
0.2 



the 

1), 
to 



Metallicity 



80k 
20k 
0.2 

Plummer 

iKroupa etall [|l993h 
range [0.1;50] M W 
IKroupa etafl [|l99ll . 

binary masses from 
100 M 
6.02 xlO 4 M 
0.36 pc 

0.53 pc 



35.8 pc 
Uniform 

Uniform in the logarithmic 
scale from 2(R 1 + R 2 ) to 50 AU 
Thermal (mod ified by 
iHurlev et all d2005l . eq. 
0.001 



1)) 



Table 4. Initial conditions for our test model 



comparison to N-body codes. For the same amount of time 
one can run multiple simulations with the MOCCA code to 
cover very wide range of initial cluster parameters. Instead 
of having one simulation from N-body, one can have hun- 
dreds of simulations from the MOCCA code and one can 
perform detail statistical analysis of the results. Addition- 
ally, MOCCA simulations give practically the same amount 
of information about the evolution of the star clusters as 
N-body codes, which makes it even more attractive. 



3 SIMULATIONS 
3.1 Initial model 

This paper is the first one of the series of papers about 
properties of BSS in star clusters. There are introduced here 
first results concerning BSS for arbitrary chosen test model. 
However, in the next papers we plan to use the MOCCA 
code to study in details how various initial conditions of star 
clusters influence on the population of BSS in general and for 
different channels of formation. We plan to run simulations 
of star clusters with various initial binary fractions, cluster 
concentrations, metallicities, different initial number of stars 
and properties of primordial binaries like eccentricities and 
semi-major axes. However, for the need of this paper, which 
is a test of the MOCCA code, we chose only one model, 
which shows the ability of the MOCCA code to follow the 
evolution of BSS population in the star cluster. 

As the test (initial) model we selected a model with 100k 
objects (80k single stars (N s ) + 20k binaries (Nb)). Initial 
conditions used in this paper are summarized in Tab. [4] 

Masses for single stars wer e chosen according to the 
initial mass function described bv lKroupa et al.ll|l993h . with 
following parameters: 



€(M) = 



M 
M~ 



-1.3 



if M < 0.5 
if M > 0.5 



(1) 



where 0.5 is so-called brake mass. Minimum star mass was 
set to O.IM© and maximum to 50Mq. All stars are assumed 
to be on the zero-age main sequence (ZAMS) in the begin- 
ning of the simulati on. Masses of binaries are chosen from 
l|Kroupa et al.lll99ll . eq. 1): 



M(X) = 0.33 



(1-X) 



0.04(1 



-(1-X) 2 



X) - 25 1.04 1 

(2) 

where X is randomly chosen from to 1. The Initial Mass 
Function (IMF) algorithm for binaries in MOCCA works in 
such a way that first it draws randomly mass for the whole 
binary according to the Eq.[5] Then, the mass ratio is drawn 
from a uniform distribution, and the whole mass is split into 
two separate masses keeping sum of the masses intact. The 
IMF generated a model with the total mass of 60200 Mq 
and with mass of the most massive binary of 89 Mq . 

Tidal radius was set to 35.8 pc. Initial observational 
core radius was set to 0.36 pc and observational half-mass 
radius to 0.53 pc. Observational core radius (r c i) is defined 
as a distance at which the central surface brightness drops 
by a half and half-light radius (ru) is defined as a radius 
within which half of the star cluster luminosity is contained. 
Metallicity was set to 0.001, which is typical for GCs metal- 
licities. 

Semi-major axes were uniformly chosen from the loga- 
rithmic scale from 2(Ri + R2) AU to 50 AU, where Ri, R2 
are stellar radii of binaries' components. Orbital periods dis- 
tribution of binaries are almost uniformly distributed, in the 
logarithmic scale, from ~ 10° up to 10 5 days. 

Eccentric ity (et) is chosen from thermal distribution 
l|Heggid Il975h but with one exception. The procedure to 
draw eccentricities modifies eccentr icity acc o rding to tidal 
evolution (equations 3b and 4 from iKroupal ()l995l l ). If the 
drawn semi-major axis is less than 5 times the ZAMS radius 
of the primary star, then it is assumed that a merger would 
happen in the beginning of the evolution, so the new set of 
parameters is altered, according to the following equation: 



Ine, = 



Rner 



+ lnej, 



(3) 



where A = 28, x ~ 0.75, and e; is the eventual eccentric- 
ity for the primordial binary. This procedure prefers small 
and very small eccentricities. Density plot with initial semi- 
major axes and eccentricities is shown in Fig. [T] It shows 
that the substantial fraction of binaries have indeed eccen- 
tricities ~ 0.0 and relatively small semi-major axes (around 
lO 5 - 1O 1O R ). 



3.2 Statistical fluctuations 

Fig.[2]shows the number of BSS from 5 simulations with the 
same initial conditions to point out how intrinsic method 
fluctuations influence the overall BSS population. BSS are 
the main subject of this paper so statistical fluctuations of 
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Mean BSS number 




2000 4000 6000 8000 10000 12000 14000 



Time [Myrs] 

Figure 1. Density plot of initial semi-major axes and eccentrici- Figure 3. la and 2a errors for the mean number of BSS from 

ties. Semi-major axes are shown in the log [R/Rq]. The substan- 5 simulations with the same initial conditions but with different 

tial fraction of binaries have eccentricities very small or equal to initial random number (seed) 

0.0 and relatively small semi-major axes (around 10 ' 5 — 10 10 Rq). 



Number of BSS from simulations with different initial seeds 

80 I 1 1 i i 1 1 

test model .seed =30 

seed =126 — 
seed =230 



70 - seed =1078 " 




1 1 1 1 1 1 ' 1 

2000 4000 6000 8000 10000 12000 14000 

Time [Myr] 

Figure 2. Total number of BSS from 5 simulations with the same 
initial conditions but with different seed values. Test model used 
in this paper was initialized with the seed=30. 



their number were used here to show the strength of the fluc- 
tuations. In Fig. [2] all plots have similar characteristics. The 
first peak of BSS is for ~ 1 Gyrs, then number of BSS drops 
and the second main peak of number of BSS lies around 
4 Gyrs. After that, it drops more or less the same for all 
simulations. 

Fig.[3]shows the mean number of BSS calculated at each 
100 Myrs from these 5 simulations together with la and 2a 
errors. Fig. and Fig. show that different seeds do not 
change characteristics of BSS population significantly. For 
almost the whole simulation la is about ± 5 BSS. In the 
worst case the standard deviation from the mean value is 
for the time ~ 5.3 Gyrs (about ± 10 BSS). 

In general all features in the number of BSS which are 
about ± 5 BSS can be the results of the statistical fluctua- 
tions. They should not be considered in the discussion when 
different models are compared. General trends in number of 
BSS are preserved. The first peak in number of BSS (~ 1 
Gyrs) and the second more extended one (~ 4 Gyrs) are 
present more or less at the same time for all 5 simulations. 
In the Fig.[3]there are shown fluctuations for all BSS from all 



channels of formation (described in Sec. I4.1.T]) . but la and 
2a fluctuations for each channel of formation separately look 
very similar. 



4 RESULTS 

In the first subsection we present a description of different 
channels of formation of BSS, together with their physical 
properties. In the next subsection there is a discussion about 
global parameters of BSS. 

4.1 Channels of formation of BSS 

In this chapter we introduce definitions of all channels of 
formation of BSS in the MOCCA code. Then, we describe 
their physical properties and discuss how BSS change their 
types during the star cluster evolution. Finally, we compare 
results of the MOCCA code with the old version of the code 
(without Fewbody) taking into account BSS population. 

All stellar evolution processes described in this paper 
are implemented only in single-star evolutio n (SSE) and 
binary-star evolution (BSE) part of the code |Hurlev et alj 
2000, 2002) and there are no additional evolutionary proce- 
dures implemented in the MOCCA code. 

4-1.1 Channels definitions 

A star is recognized as a BSS when it is a main sequence 
star and has a mass higher than the turn-of-mass by at least 
2% (to have the same condition as in J. Hurley's N-body 
simulations and t o be able to compar e our simulations with 
N-body ones, see lGiersz et alj l|201ll ) for details). Similarly 
for observations, a star which is too close to the turn-off mass 
cannot be considered as BSS because of the observational 
errors in determination of the stellar magnitudes. 

There are two main channels of BSS formation pre- 
sented in the literature. The first one is a mass transfer and 
the second one is a collisional channel of formation. We ex- 
pect that BSS created by a collision between two field stars, 
or between a binary and a field star, or between components 
of a merging binary will have different spatial distributions. 
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The details of formation processes may also vary within one 
channel. Thus, instead of having just one collisional channel 
of formation, there will be more specific channels for distinct 
scenarios. In our simulations we can trace the history of each 
object, so it is possible to find out what was the real cause 
of the creation of BSS. 

In this paper, for the simplicity of the discussion, we 
divided BSS formation channels into two general categories: 
stellar evolution and dynamical. 

There are three channels of formation of BSS which 
we include to the stellar evolution category. The first one, 
Evolutionary Merger (EM), describes a scenario when two 
stars from one binary merge into one star because of stellar 
evolution only. The second channel is called Evolutionary 
Mass Transfer (EMT). It describes a situation when there 
is some mass transfer in a binary, so that the mass of one of 
the stars overcomes the turn-of-mass. In this case the stellar 
evolution does not lead to a binary merger. This scenario is 
described in the literature as one of the two main channels 
of formation - mass transfer. The third channel is called 
Evolutionary Dissolution (ED). It describes a scenario when 
the stellar evolution leads to a disruption of a binary (e.g. 
SN explosion) with some mass accretion by the companion, 
which in consequence becomes a BSS. 

For EMT BSS it is still possible to have an evolutionary 
merger event after some time, but we did not observe such 
scenarios of changing BSS types from EMT to EM in our 
simulations (see later Sec. I4.1.5[) . Additionally, it is worth 
to mention, that these evolutionary channels of formation 
(EM, EMT, ED) are not results of dynamical interactions, 
but rather of stellar evolution. It does not necessarily mean, 
that there were no dynamical interactions for these objects 
before. It only means, that BSS creation was detected af- 
ter performing the stellar evolution step and it is assumed 
that the evolution was the direct cause of the creation of 
BSS. Field stars can influence binaries due to soft and rarely 
strong dynamical interactions. Thus, deciding whether some 
BSS is in fact purely evolutionary is possible only after a 
careful investigation of the whole history of a particular BSS. 

Channels of formation of BSS which we include to dy- 
namical category are connected strictly to dynamical in- 
teractions and are described by the following cases. The 
channel of formation called Collision Single-Single (CSS) 
describes a physical collision between two single stars. This 
is the only channel, both from evolution and dynamical cat- 
egories, which involves only two single stars. All other chan- 
nels of formation involve at least one binary. The second 
channel called Collision Binary-Single/Binary (CBS, CBB) 
describes the scenario when there is some collision between 
any two or more stars in a binary-single (CBS) or binary- 
binary interaction (CBB). 

The rest of the channels do not in fact create a new BSS 
but rather describe the change of BSS type (see Sec. 14.1.5]) . 
Exchange Binary-Single/Binary, corresponds to the situa- 
tion when BSS changes its companion in a binary, or be- 
comes a single star, or goes into a binary. Again, EXBS 
means an exchange event in a binary-single dynamical inter- 
action and EXBB means an exchange in a binary-binary in- 
teraction. The last dynamical channel of formation is called 
Dissolution Binary-Single/Binary and corresponds to the 
scenario when BSS was present in a binary, which was dis- 
rupted by a binary-single dynamical interaction (DBS) or 
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EM 


Evolutionary Merger 


EMT 


Evolutionary Mass Transfer 


ED 


Evolutionary Dissolution 


CSS 


Collision Single-Single 


CBS 


Collision Binary-Single 


CBB 


Collision Binary-Binary 


EXBS 


Exchange Binary-Single 


EXBB 


Exchange Binary-Binary 


DBS 


Dissolution Binary-Single 


DBB 


Dissolution Binary-Binary 



Table 5. Abbreviations of channels of formation of BSS and type 
changes 

binary-binary interaction (DBB). EXBS, EXBB, DBS and 
DBB cannot be initial types for BSS. Initial BSS type can 
be EM, EMT, ED, CSS, CBS or CBB, and only later BSS 
can change its type into another one. 

Similar distinction for BSS ch annels of forma tion and 
type changes are introduced also bv lLeonardl (|l996T K Abbre- 
viations of ours channels of formation and type changes one 
can find summarized in Tab. [5] 

4-1.2 Evolutionary mass transfer BSS 

For BSS from a one channel of formation there could be dif- 
ferent physical processes which are responsible for creating 
such a BSS. Thus, in this subsection BSS from the channel 
EMT will be split further into several subgroups and their 
properties will be discussed. 

Fig. [4] shows properties of BSS from the channel EMT 
divided into 3 separate subgroups: short, medium and long 
period EMT. Each subgroup corresponds to a different ini- 
tial semi-major axes range (see red, green and blue points on 
the top panels in Fig. [4]). Additionally, there are introduced 
two subgroups: dormant and strong BSS. They actually con- 
sist of points from first 3 subgroups. Dormant BSS are BSS 
for which there is some delay between the last mass transfer 
and a time when star actually exceeded the turn-off mass 
(for details about lifetimes of BSS see Sec. I4.2.3|l . The time 
when the last mass transfer occurred is called the creation 
time of a BSS, because such an event is the real cause re- 
sponsible for the creation of a BSS. The time when a BSS 
actually exceeded the turn-off mass is called the recognition 
time. Furthermore, the time range between the creation time 
and the termination time of a BSS is called the total life- 
time of a BSS. The time range between the recognition time 
and the termination time of a BSS is called the effective 
lifetime of a BSS. Strong BSS are those which had at least 
one strong dynamical interaction before BSS creation time 
(binary's properties had to be changed at least by 10%, see 
later Sec. 11X6)) . 

In the top two panels in Fig. [4] one can see how semi- 
major axes change from the initial one (at the time T = 0) 
to the semi-major axes at the recognition time. At the recog- 
nition time semi-major axes create two separate trends (the 
top-right panel in Fig. [4]). The long period EMT belong to 
the first one and short and medium period EMT to the sec- 
ond one. After a detailed analysis of the history of formation 
of these groups of EMT it came up that there are two dif- 
ferent scenarios responsible for their creation. 
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Figure 4. Properties of BSS from the channel EMT divided into 3 separate subgroups: short (red), medium (green) and long (blue) 
period EMT. Additionally, there are introduced here two subgroups: one with dormant EMT and one with EMT which had at least one 
strong dynamical interaction before becoming BSS (see text for their definitions). The top- left panel shows semi-major axes [R©] of initial 
binaries, in which EMT were created later on, with positions on the X axis set according to BSS recognition time. The top-right panel 
shows semi-major axes [Rq] of EMT in the time of their recognition. In the middle panels there are semi-major axes and eccentricities 
similarly at the time T=0 and at the time of BSS recognition. On the two bottom panels there are the mass ratios at the time T=0 
(bottom-left), and at BSS recognition time (bottom-right). On the bottom-right panel there is also one fit curve (M turn _ y j /M^rj), 
described in details in the text. In general, left panels show initial properties of binaries (T=0), where later on BSS were created, and 
right panels show properties of BSS at the recognition time. 



Long period EMT 

A donor star at some point starts to leave the main sequence, 
enters the giant branch and ejects some mass through stel- 
lar winds. A star, which later on becomes a BSS, gains 
mass of about a few ~ O.OIMq, which is not enough to be- 



come a BSS. The donor star goes further through the AGB 
phase, which means that it ejects outer layers at a rate of 
about 10~ 6 — 1CP 7 Mq per year. This stage can last at most 
few Myrs before the donor star becomes a white dwarf sur- 
rounded by a planetary nebula. The future BSS increases its 
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mass by about < O.IM© through stellar winds and by mov- 
ing through the ejected envelope. Initial eccentricities are 
typically larger than 0.1 (up to 0.9) which makes the mass 
transfer a bit easier. Although a circularization occurs, it 
does not change eccentricities to 0.0 for all binaries. Some 
of them at the recognition time are still larger than 0.0. 
The donor star becomes a white dwarf, but the star which 
gained the mass in many cases still does not exceed the turn- 
off mass. Even if the mass of the star, which gained mass, is 
close to the turn-off point, it takes some time for the BSS to 
actually exceed the turn-off mass. This is why long period 
EMT are in majority dormant BSS. It takes sometimes even 
several Gyrs until the star finally exceeds the turn-off mass 
and becomes BSS. What is also important, is that almost 
all of long period EMT are short living BSS (see Sec. 14.2. 3ft . 
It means that BSS gained only a small amount of additional 
mass from an evolving companion and exceeded the turn-off 
mass, as a MS star, just for the short period of time be- 
fore it left the MS (at most a few hundred Myrs for time 
T > 2 Gyrs). Lifetimes of BSS are described in details in 
the Sec. [4^231 

Semi-major axes of binaries with long period EMT are 
getting smaller with the recognition time (see the top-right 
panel in Fig. 0}. It is caused by the fact that the mass trans- 
fers are getting smaller with time too. In order to transfer 
enough mass to create a BSS, binaries need to have smaller 
semi-major axes. 

Many of the long period EMT are also BSS which had 
some strong dynamical interactions before becoming a BSS. 
It is most likely because of the fact that they have large semi- 
major axes and it is simply easier to change their properties 
by at least 10% in comparison to the short period EMT. 

Short and medium period EMT 

Short period EMT are created by a mass transfer in a form 
of a Roche lobe overflow. When a donor star starts to evolve 
from the main sequence its radius increases. At some point 
the binary becomes semi-detached and a conservative mass 
transfer begins. If the more massive star loses mass, then 
the semi-major axis of the binary is getting smaller. The 
mass transfer lasts typically for several hundreds of Myrs 
until the mass ratio flips in favor of the future BSS. For 
some cases, when the mass transfer continues also when the 
more massive star becomes less massive, the semi-major axis 
increases. Thus, in the top-right panel in Fig. 2] the semi- 
major axes are more scattered than initial ones. The donor 
star goes through the Hertzsprung gap phase, then the giant 
branch phase. The mass transfer stops and the donor star 
ends up as a helium core star. In many cases after some time 
such a binary merges into one star, which is not a BSS any 
longer. This mechanism is possible only for compact binaries 
(see the top-right panel in Fig. [4|. 

Semi-major axes of short period EMT at the recognition 
time decrease with time (from several hundreds Rq to about 
10 R© after about 1 Gyrs). It is caused by the fact that 
for less heavier MS stars also radii are smaller. Thus, when 
one of the stars of a binary leaves the main sequence, the 
other one has to be sufficiently close to create the Roche 
lobe. Additionally, as cluster evolves, there are less and less 
binaries which have parameters needed to create the Roche 
lobe overflow. 



On the top-right panel one can see that medium period 
EMT (green points) shrink their semi-major axes so that 
they overlap with short period EMT. In turn, short period 
EMT (red points) for time < 500 Myrs increase their semi- 
major axes. These are the only subgroups of EMT which 
change their semi-major axes that much from the begin- 
ning of the evolution. It is caused by the fact that for the 
medium period EMT, during a heavy mass loss through stel- 
lar winds, some mass is ejected outside the binary so that 
the overall binary mass is smaller. A decrease of the mass 
of the binary causes a shrinking of the semi-major axes. 
Eventually, a semi-major axis decreases to such a value that 
the short semi-detached phase in the binary starts. Mass is 
transferred via the Roche lobe overflow to the future BSS 
(this time more mass than the mass gained from previous 
stellar winds) and the star becomes a BSS. In turn for short 
period EMT, when a mass transfer occurs via the Roche 
lobe overflow in the semi-detached phase, the overall mass 
of a binary is conserved. At some point, the star which was 
loosing mass, becomes less massive (Algol-type star). The 
mass ratio inverts, but the donor star which is less massive, 
still looses its mass. Because of the conservation of the mo- 
mentum, the semi-major axis increases. As a result, medium 
and short period EMT overlap each other at the recognition 
time. 

The gap between short and long period EMT is not 
caused by the lack of binaries. There are still binaries for 
the whole range of semi-major axes, and there are still mass 
transfers through mechanisms described above. However, ei- 
ther stellar winds are not large enough to exceed the turn-off 
mass, or stars in a binary are too far of each other to create 
the Roche lobes overflows. In these cases the mass trans- 
fers are not sufficient enough to create EMT BSS. For some 
other cases the mass transfers are significant (> O.IM©) 
but the overall mass of the star does not exceed the turn-off 
mass. The turn-off point during the whole simulation does 
not decrease enough to allow these stars to become a BSS. 

Circularization of the orbits 

On the middle panels in Fig. [3] there are eccentricities 
and semi-major axes for all subgroups of EMT. Almost 
all eccentricities for short and medium period EMT were 
changed due to binary stellar evolution to = 0.0 (circular- 
ization). Circulariza tion for them is caused by tidal forces 
l|Hurlev et al. 2002). Only a few of them are not entirely 
circularized. These short period EMT which at the recog- 
nition time have eccentricities > 0.0 are created during the 
first few dozens Myrs, so there is simply not enough time to 
circularize their orbits before BSS phase. 

Eccentricities for long period EMT became smaller, 
some of them are changed to 0.0. Circ ularization for long 
period EMT is caused by stellar winds (|Hurlev et al.112003 . 
Eq. 15). For the widest binaries (> 1O 3 R0) orbits have still 
high eccentricities. They decrease in general, but not en- 
tirely. 

Masses 

Mass ratios for subgroups of EMT are shown in the two 
bottom panels in Fig. [4] These mass ratios are for all cases 
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calculated as mass of a BSS over a mass of a companion star 
(the bottom- right panel). On the bottom-left panel there are 
initial mass ratios for binaries, which later on created EMT. 

Initial mass ratios are mostly spread for short period 
subgroup of EMT (from 0.25 up to 1.0). Short period 
EMT have the mass transfers through the Roche lobe over- 
flows. When stars are close to each other, the mass transfer 
through the Roche lobes is more efficient. The companion 
looses a lot of mass, mass transfer lasts longer, and the com- 
panion cannot create a degenerated core. Thus, the mass ra- 
tios for such BSS are more spread. In the majority of cases 
the companion becomes a helium core star. Only in several 
other cases the companion evolves to WD. Thus, many of 
short period EMT at the recognition time have high mass 
ratios (up to ~ 12). 

Initial mass ratios (at the time T=0) for medium and 
long period EMT occupy similar regions: from ~ 0.5 up 
to 1.0. The mass ratios at the recognition time drops from 
about 4 (at the time T=0) to about 1.5 (at the time 
10 Gyrs). They have a narrow range at the recognition time 
because the companions evolve to WDs in almost all cases. 
Long period EMT can gain only small amount of matter 
from a companion (typically 0.1 Mq after 500 Myrs). Thus, 
for long period EMT mass ratios at the recognition time cre- 
ate a distinct trend (see the bottom-right panel in Fig. Q. 
We plot there a black curve based on the following equation: 



BSS orbital periods 



q = M P 



-off/M w 



(4) 



where M tU m-of f is the turn-off mass and M wd is a mass 
of a white dwarf (|Chernoff fc Weinberdll990l . Tab. 1). This 
equation shows how mass ratio changes through the entire 
simulation for long period EMT. The masses of BSS are 
just slightly larger than turn-off mass, thus in the nominator 
there is Mtum-off- WDs are companions in the long period 
EMT, thus in the d enominator there are masse s of WDs 
calculated based on IChernoff fc Weinberg! <|l990l . Tab. 1). 
One can see that the black curve in the bottom-right panel 
in Fig. |4] reproduces long period EMT mass ratios very well. 
They have a narrow range and predictable values through 
the entire simulation. 



4-1.3 Evolutionary merger BSS 

Fig.[5]shows last orbital periods before a star was recognized 
as a BSS as a function of the recognition time. EM BSS 
are compact objects, with orbital periods ~ 1 day or less 
and circular, or almost circular, orbits. They were created 
also from the compact binaries with semi-major axes, at the 
T = 0, from about 10 5 up to about 1O 12 R (see bluest 
points in the density plot in Fig. [1} and with eccentricities 
~ 0.0. These are expected properties for EM BSS, because 
in order to have an evolutionary merger, there has to be a 
very close binary. 

For the whole simulation there were observed only two 
exceptions for EM channel with the orbital periods ~ 10 4 
days. After checking the history of these long period EM 
it came up that these were binaries which were changed by 
interactions, got very high eccentricities (> 0.99) and after 
some time the binary merged into one star creating BSS. 

Fig. [6] shows masses of BSS at the recognition time for 
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Figure 5. BSS orbital periods at the recognition time for different 
channels of formation 
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Figure 6. Histogram with masses of BSS divided by the turn-off 
mass at the recognition time for different channels of formation 



3 major channels of formation: EMT, EM and CBS+CBB. 
The majority of EM BSS have masses which are just 
slightly above the turn-off mass. Almost all of such EM 
(1 — 1.1 Mtum-off) are the dormant EM (see Sec. 14.2.3)) . 
When dormant EM are detected, they have masses just 
slightly larger than the turn-off mass. For larger EM masses 
(> 1.1 Mtum-off) there are less and less EM. In order to 
create an EM with mass close to 2 Mtum-off both com- 
ponents of a binary would have to have masses close to the 
Mtum-off - There is simply less compact binaries with initial 
mass ratios close to 1. For compact binaries IMF generates 
rather different masses of components. However, it is still 
possible to have EM BSS with higher mass ratios, which 
eventually create heavier BSS (up to 2 times the turn-off 
mass). 

In Fig. [5] one can see that EM BSS consist actually 
of two subgroups. The first group represents EM with the 
recognition times ^ 3 Gyrs for which orbital periods clearly 
decrease with time. At the beginning of the simulation or- 
bital periods of binaries, which later on become BSS, are 
slightly larger than 1 day. Just before 3 Gyrs orbital peri- 
ods, which give rise to the number of BSS, shorten to about 
10 -0 ' 5 days. Around 3 Gyrs they increase and stay around 
1 day for the rest of the star cluster evolution. We carefully 
checked the history of the creation for many BSS of these 
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two groups and we were able to find two different physical 
scenarios responsible for the creation of BSS. 

In the first group of EM 3 Gyrs) the typical physical 
process responsible for the creation of EM BSS is a merger in 
a binary due to the Roche lobe overflow. Typical formation 
scenario starts with the binary with masses like e.g. mi = 
3.3M0, 77i2 = O.5M0, the semi-major axis a = 6.2R©, and 
the eccentricity e = 0.004. The semi-major axis decreases 
with time but only slightly, the eccentricity circularizes to 
0.0, and after a few hundreds Myrs, when the heavier star 
increases its radius, a semi-detached phase starts. It lasts 
less than 1 Myr and finally the binary merges into one BSS 
star. 

The typical scenario for EM BSS created after 3 Gyrs 
involves the magnetic braking which eventually leads to a 
merger event in a binary. The magnetic braking affects stars 
with m asses less than 1.25M0 (for details see iHurlev et al.l 
(2002)). Indeed, at the time T ~ 3 Gyrs the turn-off mass 
equals 1.25M0 and the magnetic braking starts to work for 
both components in binaries with MS stars. Typical for- 
mation scenario for EM BSS recognized after 3 Gyrs is as 
follows. The initial binary properties are e.g. mi = O.88M0, 
m 2 = O.13M , a = 8.89R©, e = 0.003. The masses and 
the semi-major axis are in general smaller in comparison 
to the binaries which created EM through the Roche lobe 
overflow. The semi-major axis decreases after 8 Gyrs up to 
a — 3.41R©, the eccentricity circularizes to 0.0. The stars 
merge into one object, without detecting the semi-detached 
phase by the MOCCA code. This and larger decrease of the 
semi-major axis (through the magnetic braking) are two dif- 
ferences between these two groups of EM BSS. 

Additionally, these two groups of EM BSS still over- 
lap because magnetic braking works for binaries with stars' 
masses below 1.25Mq, and there are such binaries in the 
system before 3 Gyrs. It is more complicated also because 
of the existence of a dormant EM BSS for which the merger 
event can be significantly earlier than the recognition time 
(see Sec. 14.2. 30 . EM BSS created by the Roche lobe overflow 
are created mainly for the more massive stars in the first few 
Gyrs and those created by the magnetic braking are created 
mainly for the later stages of the star cluster evolution. But 
still, there are some EM created by magnetic braking before 
3 Gyrs (for masses less than 1.25Mq), and there are some 
Roche lobe overflow EM BSS created after 3 Gyrs. It will 
be interesting to see if a similar division is present also for 
different initial binary properties. We plan to check this in 
details in the next paper. 

4-1-4 Collisional BSS 

Collisional BSS (CBS and CBB) are created in strong dy- 
namical interactions between binaries and single stars. Fig. [5] 
shows orbital periods of binaries at the recognition time. 
Before core collapse (< 6 Gyrs) CBS and CBB were cre- 
ated from binaries with rather large orbital periods. It is 
caused by the fact that before the core collapse, when den- 
sity is smaller, only very wide binaries have probabilities 
large enough to have dynamical interactions. Additionally, 
some of such primordial binaries are already in the core due 
to the initial conditions. During and after the core collapse, 
when the core is getting denser, there are more CBS and 
CBB in general. They are created also from more compact 



binaries with orbital periods of about a few up to 10 days. 
After the core collapse, the density in the core is high enough 
to create BSS, even from close binaries. We checked also the 
initial properties of the binaries, which later on created CBS 
and CBB, and found no correlations with the orbital peri- 
ods. Properties of such primordial binaries are uniformly 
distributed for the whole range of the semi-major axes and 
eccentricities in the histogram in Fig. [1] 

Fig-E]shows masses of BSS divided by the turn-off mass. 
CBS and CBB have roughly uniformly distributed masses 
from 1 up to 2 times the turn-off mass. There is no dis- 
tinct peak in the distribution of the masses for this channel. 
It means that CBS and CBB BSS are formed from various 
MS stars and there is no particular preference which of them 
collide in dynamical interactions. The only limitation is that 
two MS stars can form BSS with the mass equal to at most 
twice the turn-off mass. However, there are still several CBS 
or CBB with masses exceeding the turn-off mass more than 
two times. These BSS were created by rather complex sce- 
narios. In some cases there were more than one merger for a 
star which became a BSS. In other scenarios, there was an 
additional mass transfer from the companion to the already 
existing BSS, which gave BSS enough mass to exceed the 
turn-off mass more than two times. 

Initial binaries' mass ratios for CBS and CBB channels 
are uniformly distributed up to 2.5 for various eccentricities, 
up to 0.8. CBS and CBB become important after several 
Gyrs, and during that time there are actually no more bi- 
naries with mass ratios larger than about 3. CBS and CBB 
binaries' mass ratios more or less cover overall mass ratios 
of all binaries in the system at a given time. There is no 
visible preference which binaries could create CBS or CBB 
in collisions during dynamical interactions. 

Apart from to the fact that CBS and CBB were created 
in dynamical interactions, there are no distinct physical pro- 
cesses or typical formation scenarios for creating these types 
of BSS. Collisions occurred for binaries with various semi- 
major axes and eccentricities. Some dynamical interactions 
were simple, other ones were resonant. For some other cases, 
during the dynamical interactions, the collision occurred for 
a binary after it already passed the other object (a binary or 
a single star) at the pericenter. The binary, going through 
the pericenter, gained a very high eccentricity and a collision 
in a binary occurred while it was already moving away from 
the other object. But still such dynamical interactions were 
not resonant because the scattering objects run through the 
pericenter only once. 

4-1.5 Types changes 

Fig. shows the number of BSS as a function of time. One 
can see how different channels of formation of BSS change 
their significance during the evolution of the star cluster. 
EMT channel is the most active in the beginning of the sim- 
ulation. In the beginning, there are simply more binaries 
which are able to start a mass transfer. First EMT BSS are 
formed in very massive binaries, which means that these bi- 
naries evolve fast from MS, there occur violent mass trans- 
fers which create short living BSS (see Sec. I4.2.3|) . Stellar 
evolution for massive stars is very fast and after several hun- 
dreds of Myrs, when EMT gains its peak, the EMT channel 
becomes less and less significant as the cluster evolves. The 
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number of BSS created by the EMT channel continuously 
decreases but stays active for the whole simulation. The 
number of EMT BSS decreases because it depends mainly 
on the initial properties of the binaries. Depending on the 
initial mass, a star leaves MS at a different time. Thus, dif- 
ferent binaries starts the EMT BSS phase at different times. 
As the evolution of the star cluster continues, the number 
of EMT BSS continuously drops, because there are less and 
less binaries which can start mass transfers and create EMT 
BSS. 

The number of EM BSS in Fig. [7] starts to increase af- 
ter 1 Gyrs, when a violent evolution of most massive stars is 
over. A binary, depending on the physical process responsi- 
ble for the creation of BSS (the Roche lobe overflow or mag- 
netic braking), needs some time to actually merge. Thus, 
the first increase of EM BSS appears after around 2 Gyrs. 
The peak of the number of EM falls to around 4 Gyrs and 
after that it decreases. The peak of the number of EM BSS 
is most likely caused by the fact that both scenarios for 
EM creation, the Roche lobe overflow and magnetic brak- 
ing, work at that time most efficiently. The magnetic brak- 
ing (see Sec. I4.1.3P works for stars with masses smaller than 
1.25Mq . Thus, around the time T = 3 Gyrs it starts to work 
for both components in binaries with MS stars. Because of 
that, it is possible to create mergers more efficiently and also 
from slightly wider binaries. Thus, there is a distinct peak in 
the orbital period values for EM around 3 Gyrs (see Fig. [5} 
and an increase of the number of EM around 3-4 Gyrs (see 
Fig. [7J . There are simply more binaries with proper initial 
parameters to create EM BSS through the magnetic brak- 
ing. 

CBS and CBB are created in star cluster from the be- 
ginning, but significantly more of them are created during 
a core collapse and after that (> 6 Gyrs). CBS and CBB 
are created mainly in the core up to the half-mass radius 
(see Fig. I12|l . at the time of the core collapse (8 Gyrs), and 
afterwards, when the core is dense, they become one of the 
most active channels of formation of BSS. 

In Fig. [7] one can see that there is only one BSS, which 
was created by the CSS channel (about 3 Gyrs). Just to 
remind, the CSS channel is the only channel which creates 
BSS from a collision between two single stars. All the other 
channels involve one or more binaries. It means, that actu- 
ally almost 100% of BSS, at least for this test simulation, 
were created from binaries. Even if the CSS channel is not 
important for our test model, it does not necessarily mean 
that it will be not important for other simulations with dif- 
ferent initial conditions. Perhaps for star clusters with higher 
concentrations there will be more CSS? 

We already mentioned in the Sec. 14.1.11 that types of 
BSS at the creation time can change. For example, a binary 
with BSS could be disrupted, or could exchange with some 
other star in a dynamical interaction. In such a case the 
actual type of BSS at a given time of the simulation can 
be different. Fig. [8] shows how types of BSS change during 
the simulation for different channels of formation. Just by 
quick comparison between two figures Fig. [7] and Fig. [8] one 
can see that some BSS changed their types during the sim- 
ulation, but it concerns mainly BSS created in interactions 
(CBS and CBB). In turn EMT and EM BSS did not change 
significantly during the simulation. The plots with BSS at 
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Figure 7. BSS types at the creation time as a function 
of the time 
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Figure 8. Actual types of BSS as a function of the 
time 

the creation time look like these with actual types of BSS 
for these channels (see Fig. [Jjand Fig. [HJ. 

Fig. [9] shows in details how many BSS changed their 
types. On the X axis there are types of BSS at the creation 
time and on the Y axis there is the final (last known) type of 
the BSS. One can see that only a few BSS which were formed 
as EM or EMT changed their types. There was only one EM 
which got into a binary (EM — > EXBB). There were also 
two cases where EM first, as s single star, got into a binary 
during some dynamical interaction, and afterwards such a 
binary was disrupted (EM — > DBB). There was only one case 
when EMT was disrupted (EMT -> DBB). Despite these 
interesting scenarios, from statistical point of view, BSS type 
changes for channels EM and EMT are not significant. In the 
next paper, we plan to study BSS population depending on 
many different initial conditions, both for binaries and global 
star cluster properties. It will be interesting to see how type 
changes will look like for other models. For example, for star 
clusters with higher concentrations, we expect to have more 
BSS type changes for EM and EMT channels. 

BSS from EM and EMT channels do not change their 
types significantly during the star cluster evolution. Some 
of them were created in central regions, where probability 
of dynamical interactions with other objects is higher (see 
Fig. I12[) . However, even with higher dynamical interaction 
probabilities, there are only several EM and EMT which 
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Figure 9. BSS type changes between different BSS types. On the 
X axis there is the BSS type at the creation time, and on the Y 
axis there is the final (last known) type of the same BSS. Each 
number shows how many BSS changed their types from one type 
to another. 



changed their types (see Fig. [9}. EMT do not change types 
probably because of the fact, that most of them are very 
hard binaries (see Fig.[5|, so there is needed another massive 
binary or a star and small impact parameter to disrupt or 
exchange BSS. Most likely, these binaries are hard enough 
to survive such interactions. For other EMT BSS, which 
have larger orbital periods (see Fig. [5j, the reason why they 
were not disrupted, is that these EMT live shortly (max 
~ 200 Myrs, see Sec. 14.2.3]) . so there is simply not enough 
time to have strong dynamical interactions. Some of the long 
period EMT are found also around half-mass radius and be- 
yond, where probabilities of the dynamical interactions are 
extremely small. Additionally, after the core collapse, when 
the dynamical interactions become important, there is sim- 
ply small overall number of EMT. EM do not change types 
most likely because they are single stars, and the probabili- 
ties of the dynamical interactions for them is lower than for 
EMT. Additionally, EM are not the most massive objects 
in the star cluster. Their masses are, in majority of cases, 
less than 1.2 of the Mtum-of f for time > 8 Gyrs, when the 
dynamical interactions are most important. Thus, the dy- 
namical interactions occur rather for massive WDs and NS. 
Moreover, the number of EM is small in comparison to the 
number of WDs or NS at the later stages of the star cluster 
evolution. 

The only channels which were significantly changed dur- 
ing simulations are CBS and CBB. One can see that BSS ini- 
tially created in these channels change their types to EXBS, 
EXBB, DBS, DBB. It shows that the dynamical interactions 
(exchanges and dissolutions) starts to play a significant role 
for CBS and CBB in the later stages of the star cluster evo- 
lution (see Fig.|SJ). It correlates with the core collapse event, 
which starts at about 6 Gyrs, and after the core collapse, 
at ~ 8 Gyrs. Most of CBS and CBB, after the dynamical 
interactions which created them, are still inside binaries. For 
the binaries there is a higher probability that dynamical in- 
teraction will occur. Additionally, BSS in binaries are heavy 
objects and thus, they sink to the center of the star cluster, 
where dynamical interactions are even more frequent. 14 of 
the total 36 CBS changed their types, which gives ~ 39% 
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Figure 10. Histogram showing the number of BSS, binned for 
1 Gyr, which might be created by induced mass transfer or in- 
duced evolutionary merger (see text for explanation). Each BSS 
is taken into account only once, according to its recognition time. 
The value is divided by the number of all EM or EMT BSS in the 
particular time of the simulation, to find out whether this pro- 
cess is significant. Note: if there are no histogram bars for some 
time ranges, it means that there arc no BSS with possible induced 
mass transfer (or merger) created for these time ranges, but there 
are still some BSS in the system in general. 



efficiency. Also, 18 of the total 58 CBB changed their types, 
which gives ~ 31% efficiency. 

During the entire simulation there were created overall 
476 BSS. The most active channels of formation were EMT 
and EM. There were created 231 EMT BSS (49%) and 149 
EM BSS (31%). There were 95 BSS from the CBS and CBB 
channels (20%), and there was just one BSS created by the 
CSS channel in our test simulation. Despite the small over- 
all number of CBS+CBB in comparison to the evolutionary 
BSS (EM, EMT), they can be dominant channels of for- 
mation of BSS for clusters which are in the post collapse 
phase. Furthermore, for different initial conditions, we ex- 
pect to have also different number of EM and EMT. These 
two channels are a direct consequence of binary initial con- 
ditions, thus for the simulations with less initial number of 
compact binaries (see Sec. I3.1J1 we expect to have much less 
EM and EMT BSS in general. EM particularly are formed 
only from primordial compact binaries with e ~ 0.0. 



4-1.6 Possible induced mass transfer or merger 

Fig. QJj] shows EM and EMT created by a possible induced 
mass transfer or merger. Each histogram was normalized to 
1.0 by dividing BSS number by the overall EM or EMT 
number of BSS in the particular time of the simulation. 
The possible induced EM or EMT represents BSS which 
had some strong dynamical interactions involving binaries, 
before the creation time. One can call these scenarios as 
possible induced mass transfers (for EMT) or possible in- 
duced mergers (for EM). Strong dynamical interactions in 
this context is a binary-single or a binary-binary dynamical 
interaction, which changed the semi-major axis or eccentric- 
ity by at least 10%. One can see how the possible induced 
mass transfer could change its significance during the star 
cluster evolution. Induced EMT BSS in the beginning of the 
simulation are not significant because less than 5% of them 
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had some strong dynamical interactions before. Their sig- 
nificance increases with time as the core density increases. 
Later on, possible induced mass transfer increases consider- 
ably after the core collapse at 8 Gyrs, becomes important 
and concerns more than 50% of all EMT BSS. BSS from the 
EM channel increases with time too but in this case there 
is less EM than EMT. EM are created from more compact 
binaries, so there is a smaller probability of dynamical in- 
teractions for them in comparison to EMT. 

In the Sec. 14. 1.71 we present a simulation with the same 
initial conditions as for the test simulation, but with the old 
version of the code (without Fewbody). The overall number 
of EM and EMT for both simulations, with and without 
Fewbody, are essentially the same (see Fig. 1 1 1 p . It turned 
out that the number of possible induced EM and EMT for 
the simulation without Fewbody increase with time, just like 
for our test simulation. The only difference is that for simu- 
lation without Fewbody there is in general less EM and EMT 
with possible induced mass transfer. Next, we performed a 
population synthesis for the same set of primordial binaries 
as for our test simulation, to check how many of EM and 
EMT will be created (see Fig. [TTj) . We found that their num- 
ber is essentially the same as for the test simulation and the 
simulation without Fewbody. The only significant difference 
in the number of EM is around 3-5 Gyrs. We were unable to 
determine what is the cause of this difference. The number of 
EM and EMT seems to not depend on the dynamical inter- 
actions. Perhaps the induced mass transfer is not important 
after all. 

We study in this paper only one test model. A number 
of possible induced mass transfer BSS appears to be impor- 
tant only at the later stages of the star cluster evolution - 
after the core collapse. On the other hand, the numbers of 
EM and EMT are essentially the same for all three cases: 
with and without Fewbody, and for the population synthe- 
sis (so without any dynamical interactions) . It favors rather 
the conclusion that the dynamical interactions do not have 
significant influence on the creation of EM or EMT. Per- 
haps possible induced mass transfer will be more signifi- 
cant for some specific initial conditions (star clusters with 
higher concentrations?). Additionally, we checked here the 
dynamical interactions which changed binaries' properties 
at least by 10% (arbitrary chosen value). Maybe this value 
is too large and changes of parameters of binaries, which are 
smaller but frequent, should be t aken into account a s well? 
In the next paper in the series (jGiersz et alj 1201 ll ) there 
are compared MOCCA simulations with N-body ones. For 
a simulation with higher r pmax in the MOCCA code (higher 
impact parameters), there are significantly more fly- by in- 
teractions for binaries. In such a case there were also created 
many more BSS. This, in turn, suggests that the dynami- 
cal interactions could have some influence on the creation 
of BSS. For details see iGiersz et~ai1 <|201ll . Fig. 7). We will 
try to determine whether the dynamical interactions have 
indeed some implications for the population of EM or EMT 
in our future work. 

4-1.7 Comparisons with old version of the code 

Fig. [TT] shows the number of BSS from different channels of 
formation for the test simulation and for the simulation with 
the same initial conditions but with the old version of the 
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Figure 11. Comparison of the numbers of BSS for the test sim- 
ulation with the old version of the code (without the Fewbody), 
and with the simple population synthesis (without any dynamical 
interactions) . 



code (without the Fewbody). In the old version, instead of 
the Fewbody, there are used cross-sections to determine the 
outcomes for the dynamical interactions. Fig. [TT] shows also 
the number of EM and EMT obtained from simple popula- 
tion synthesis (without any dynamical interactions). 

It was mentioned before that EM and EMT BSS de- 
pend strongly on the initial binary properties. Thus, overall 
population of EM and EMT with and without the Fewbody 
is similar (see Fig. [TT}. The peak value of EM is equal to 
about 40 BSS for both simulations, and about 25 BSS for 
EMT channel. In turn, for BSS created in the dynamical 
interactions there is a huge difference. The old version of 
the code did not create even one CBS or CBB BSS because 
interactions which create mergers in dynamical interactions 
are rather complex. The old version of the code could not 
deal with them in the N-body manner. With the Fewbody 
the number of CBS and CBB outnumbers the EMT channel 
after 6 Gyrs. After 9 Gyrs it becomes also as efficient as EM 
channel. This shows how important it was to incorporate 
the Fewbody into the MOCCA code in order to follow the 
evolution of peculiar objects, like BSS. For other models of 
star clusters, for which the dynamical interactions are even 
more important, one can expect even larger differences in 
favor of the simulations with the Fewbody. 

4.2 BSS global parameters 

In this section we discuss global BSS parameters, like their 
positions in the star cluster, lifetimes, and what is the frac- 
tion of BSS in binaries in comparison to single BSS. 

4-2.1 BSS initial and final positions 

Fig. Q2] shows the initial and final positions from the cen- 
ter of the star cluster as a function of BSS recognition time. 
The initial positions refers to the distances for which the last 
merger (for all channels except EMT), or last mass transfer 
(for EMT) occurred. Additionally, in Fig. [12] to better vi- 
sualize distances, t here are shown the core r adius (r c , calcu- 
lated according to ICasertano fc Hutl (|l985l )). the half-mass 
radius (r^), and the tidal radius. 
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Figure 12. In the left panel there are BSS distances [pc] from the center, when last merger or last mass transfer occurred, as a function 
of the recognition time. In the right panel there are last known positions of BSS for different initial channels as a function of BSS 
recognition time. The core radius (r c ), the half-mass radius (r h ), and the tidal radius are also shown. 



EMT BSS initial positions extends from deep inside the 
core up to far outside of the half- mass radius (r^). Final 
positions for EMT are more or less the same. There is only 
a few EMT for which final orbits were a bit more extended. 
Only for one EMT the final position was located outside the 
tidal radius. EM BSS were created, in general, at the same 
distances as EMT. Final positions of EM are more or less 
the same as the initial ones. Only some of EM moved deeper 
to the core and some other moved near the tidal radius. 

Initial positions of CBS and CBB correlate nicely with 
the core radius (r c ). CBS and CBB before the core collapse 
were created mainly inside the core radius in the dynamical 
interactions with wide binaries. During the core collapse and 
afterwards, they were created mainly in the core, but some 
of them also outside the core up to the half-mass radius. CBS 
and CBB are very often inside binaries, so they are heavier 
than EM or EMT and sink to the center of the cluster due to 
mass segregation. Some of them, because of the dynamical 
interactions, are ejected to the halo of the star cluster. Thus, 
their final positions is far outside the half-mass radius. 

4-2.2 Bimodal distribution 

Bimodal d istribution of BSS wa s found in ma n y star clus- 
ters (|Ferraro et all (|l99l 1 19971 ); IZaggia et all (|l997r i, and 
other references in Sec. [T]). Projected distances of BSS for 
our test simulation are shown in the first panel in Fig. [13] 
for the time 4 Gyrs, and in the second panel in Fig. [T3] 
for the time 10 Gyrs. Because the number of BSS is rather 
small, in order to decrease the statistical fluctuations, the 
projected distance was calculated 10 times for each BSS. In 
Fig. [13] there are given mean numbers of BSS depending on 
the distance from the center. On the bottom label of the X- 
axis there is the distance relative to the core radius, and on 
the top label there is the distance relative to the half-light 
radius. Both, the core radius and the half-light radius, are 
observational ones (see definitions in Sec. 13. 1|) . 

In the first panel of Fig. [13] (time 4 Gyrs) one can see 
that the bimodality of BSS distribution is not present. The 
majority of BSS is located near 2 half-light radii. The BSS 
distribution continuously drops with the distance from the 
center of the star cluster. It seems that a weak bimodality 



becomes visible after the core collapse (the second panel of 
Fig. I13|) . One can see that the first high peak of the number 
of BSS is about 5 r c ; (0.5 r^i), a well defined minimum about 
18 r c i (1.8 rhi), and the second p eak in the number of BSS i s 
in the outskirts about 2.8 r M . In lFerraro et alj |l993l . I1997V I 
the plot with the bimodal distribution for M3 has wider bins 
in the histogram for the second peak value. However, in our 
plot the width of bins is the same for all distances, thus in 
our model the second peak is less visible than in these cited 
papers. 

In our next paper we plan to investigate the bimodal- 
ity for various star cluster properties and we will check for 
which star clusters and when the bimodality starts to ap- 
pear. Perhaps we will find that the bimodality is present 
only in star clusters after the core collapse? It would be a 
great tool to probe dynamical stages of stars clusters. 

4.2.3 BSS lifetimes 

Fig. [14] shows the effective (left panel) and the total (right 
panel) lifetimes of BSS for different channels of formation. 
Just to remind, the effective lifetime is calculated as a differ- 
ence between the recognition time and the termination time. 
The total lifetime is calculated from the creation time and 
the termination time. In Fig. Q3]there are also presented two 
boundary lines: the upper and the lower limits. The upper 
limit is the expected lifetime of a star on the main sequence 
with one turn-off mass, and the lower limit is the expected 
lifetime of a star on the main sequence with two times the 
turn-off mass. These boundaries show what are the theoret- 
ical minimum and maximum lifetimes for a BSS which had 
a one merger event before becoming a BSS. One can see 
that in the first panel of Fig. 1141 there are many BSS, from 
different channels of formation, for which lifetimes are signif- 
icantly smaller than the lower limit. However, in the second 
panel of the Fig. [14] there are only several such cases which 
have still smaller lifetimes than the lower limit (details are 
given later in this section). 

All BSS created in the first 1 Gyr, as expected, are 
short living BSS because these are the heaviest stars. In 
the beginning of the star cluster evolution there are many 
massive stars in binaries which drive the EMT channel to 
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Histogram of BSS projected positions for time 4 Gyrs 

'I'm 




Histogram of BSS projected positions for time 10 Gyrs 
r/r h . 




Figure 13. Histograms with BSS projected distances divided by r c i at the time 4 Gyrs, and time 10 Gyrs, for all BSS. Bimodality in 
the BSS distribution is not visible yet for the time 4 Gyrs, but a weak bimodality is visible for the time 10 Gyrs. Both radii, r c ; and r^i, 
are observational radii. 
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Figure 14. Left: the effective lifetimes of BSS [Gyrs] at the recognition time [Gyrs] for different channels of formation, calculated from 
the time range between the recognition time and the termination time. Right: the total lifetimes of BSS [Gyrs] at the recognition time 
[Gyrs] for different channels of formation, calculated from the creation time and the termination time. These two panels show that for 
many BSS there is a significant delay for a star before it can be recognized as a BSS (i.e. it exceeds the turn-off mass). The upper limit 
is the expected lifetime of a main sequence star with the mass equal to the turn-off mass, and the lower limit is the expected lifetime of 
a main sequence star with the mass equal to twice the turn-off mass. 



produce BSS through mass transfers. After a few hundreds 
Myrs, when the most massive stars leave the main sequence, 
the number of EMT starts to decrease and BSS are created 
more efficiently in other channels. Lifetimes of BSS produced 
after 1 Gyr are significantly longer than for EMT created 
earlier. The reason is that BSS created after the first Gyr 
are not as massive as in the beginning of the simulation. 
Their evolution is much slower and the BSS phase can last 
much longer. 

Fig. 1 141 shows how the BSS phase can be delayed before 
the actual detection by observations. The BSS last merger or 
the last mass transfer can happen even several Gyrs before 
the turn-off mass will decrease enough (see Fig. I15|l . so that 
star can be recognized as a BSS. Total lifetimes take into 
account delayed detection. In this case many of BSS moved 
with their lifetimes between the upper and lower limits. Ad- 
ditionally, one can see in the second panel of Fig. 1141 that 
many of BSS with corrected lifetimes also moved closer to 
the upper limit (red line). 

In the second panel of Fig. 1 141 there are about 30 more 



EMT BSS (mainly for T < 1 Gyrs) which have still total life- 
times significantly shorter than the lower limit, even taking 
into account a delayed detection. It was unexpected to have 
BSS created at later stages of the star cluster evolution with 
such short lifetimes (only ~ 10 Myrs). Some of these stars 
exceeded the turn-off mass by just a little and started a very 
short BSS phase. After just few Myrs (up to 100 Myrs) they 
evolved from MS to the red giant phase, ending the BSS 
phase. In a few other cases, instead of the red giant phase, 
there was a merger between the binary components which 
terminated the BSS phase (the resulting star was not a MS 
star). There are also two CBS/CBB with very short total 
lifetimes. For one of them, the BSS left MS very fast. For 
the second one, there occurred another merger event while 
the star was still BSS, but after that the star was not BSS 
any longer. 

Fig. [15] shows a relation between the recognition time 
of BSS and the creation time. Every point which is below 
the imaginary straight line x = y denotes BSS for which 
a merger or a mass transfer event occurred before the star 
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Figure 15. The recognition time of BSS [Gyr] vs. the creation 
time [Gyr] for different channels of formation. The plot shows 
what is the difference in time between an event which really cre- 
ated BSS and the actual recognition of a BSS. 



was detected as a BSS. One can see that many BSS from 
EM and EMT channels, and some BSS from CBS and CBB 
channels had mergers even several Gyrs before they became 
BSS. Mergers or mass transfers occurred for these stars ear- 
lier, but masses were too small, so they did not exceeded 
the turn-off mass and they had to wait for a detection as 
a BSS. This effect was not expected in common scenarios 
for creation BSS, and it was rather assumed, that mergers 
between stars create BSS immediately. This plot shows that 
delayed detection of dormant BSS is significant. For all 149 
EM BSS there were 45 dormant EM BSS (30%), for all 231 
EMT BSS there were 60 dormant EMT (26%) and for all 95 
CBS and CBB there were only 7 dormant CBS, CBB BSS 
(7%). For the total 476 BSS there were overall 112 dormant 
BSS which is 24%. 

Fig. 1 161 shows the escape time of BSS depending on the 
recognition time. It shows the escape times for the main 
channels of formation: EM, EMT and CBS+CBB. In the 
figure there are only BSS which escaped from star cluster as 
BSS. All other BSS which first stopped being BSS, and then 
escaped from the star cluster, are not taken into account in 
this plot. Also, escaped BSS in Fig. [T6]are divided into two 
groups, according to different reasons of the escape: due to 
energetic dynamical interaction (empty squares in the plot) 
and due to the relaxation process (stars in the plot). One 
can see that the number of BSS escapers increases during 
later stages of the cluster evolution, when the core starts to 
collapse (~ 6 Gyrs). The number of escapers increases most 
significantly after the core collapse (> 8 Gyrs). Additionally, 
there are significantly more BSS escapers from CBS+CBB 
channels, for which the time between the creation and the 
escape time is small (fast escapers). This is caused by the 
fact that BSS created due to dynamical interactions gain sig- 
nificant kinetic energy, which allows them to leave the star 
cluster faster. All BSS from CBS and CBB channels, except 
just one, escaped from the system due to the dynamical in- 
teractions (blue squares). In contrast, for EM and EMT BSS 
the process responsible for the escape from the system is the 
relaxation. There is just one EM case and one EMT case for 
which the dynamical interactions were the reasons of the 
escape. Relaxation do not give BSS often such a large ad- 
ditional kinetic energy as the dynamical interactions. Thus, 
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Figure 16. The recognition time [Gyr] vs. the escape times of 
BSS [Gyr] for different channels of formation. In the plot there are 
only BSS which escaped as BSS. For each channel of formation 
there are two types of points showing the cause of the escape: 
relaxation or dynamical interaction. 
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Figure 17. The ratio between the number of BSS as single stars 
and the number of BSS inside binaries as a function of time 



there are in general less EM and EMT which can leave the 
system and they need significantly more time for the escape 
(slow escapers, see Fig. I16[) . 

From the overall 476 BSS created in the simulation, 
62 escaped as BSS (~ 13%). From the total 231 EMT BSS 
created during the whole simulation, only 8 of them escaped 
as BSS (~ 3%), and from the total 149 EM BSS, only 14 
escaped as BSS (~ 9%). However, BSS escapers are more 
important for CBS and CBB channels because for the total 
95 of them, 40 escaped as BSS (~ 43%). CBS and CBB 
escapers become even more important after the core collapse 
(> 8 Gyrs), when there were created 50 CBS and CBB and 
30 of them escaped as BSS (60%). For our test model on 
average over 10% of BSS escape from the cluster, so it seems 
that escape BSS process could be important. Thus, one can 
try to search for them in tidal tails of post-collapse stars 
clusters (especially slow escapers). 



4-2.4 BSS in binaries 

The ratio between the number of single BSS and the number 
of BSS in binaries as a function of time is shown in the 
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Fig. 1171 In other words, plot shows how many times there is 
more BSS as single stars than BSS in binaries. 

The ratio starts favoring BSS in binaries. The main BSS 
channel responsible for that is EMT. The mass transfer cre- 
ates BSS in primordial, massive binaries, but the EMT chan- 
nel dominates only for the first ~ 2 Gyrs of the star cluster 
evolution. After that time EM start to play a bigger role, 
as single BSS. The number of BSS from EM channel contin- 
uously increases, while EMT BSS continuously decreases, 
until there is 2-3 times more BSS as single stars than in 
binaries. The shape of the plot, especially its peaks, corre- 
sponds nicely to the peak of the EM channel (see Fig. [7J). 
When EM channel starts to drop after ~ 4 Gyrs, BSS in 
singles to BSS in binaries ratio starts to drop as well. The 
ratio drops to ~ 1.5. It is interesting that it stays around 
this value for the next several Gyrs. Just at the end of the 
star cluster evolution (> 12 Gyrs) the ratio flips again, and 
there are more BSS in binaries again. This is caused by the 
fact that at the end of the star cluster evolution there are 
more binaries left, because of their heavier mass. 



5 DISCUSSION AND SUMMARY 

We described here an improved code for simulations of star 
clusters, called MOCCA, which stands for MOnte Carlo 
Cluster simulAtor. It combines two very distinctive ap- 
proaches for the star clusters simulations: the Monte Carlo 
method and the Fewbody - a direct integration method. The 
Fewbody is used to perform small N-body scattering experi- 
ments between binaries and single stars and between two bi- 
naries. This gave the current version of the code very needed 
features. From now on, one can follow the evolution of exotic 
objects formed in complicated dynamical interactions, like 
BSS. The MOCCA code was used here to investigate a popu- 
lation of different channels of formation of BSS in star cluster 
environment with 100k initial objects (80k single stars and 
20k primordial binaries). This model was chosen as a test 
model to check the ability of the code to follow the evolu- 
tion of BSS. BSS were divided in this paper into two major 
categories: evolutionary and dynamical. To the evolution- 
ary category of BSS we include evolutionary mergers (EM), 
evolutionary mass transfers (EMT) and evolutionary disso- 
lutions (ED). EM is created when two stars from a binary 
merge as a result of the binary stellar evolution (without dy- 
namical interaction with other stars in the system). EMT is 
a BSS which is created due to a mass transfer, which allows 
the star to exceed the turn-off mass. ED is a BSS which is 
created when a binary dissolves during the supernova explo- 
sion. To the dynamical category of BSS we include all BSS 
which were created as a result of the dynamical interactions 
between two single stars (CSS), or between single star and 
binary (CBS), or between two binaries (CBB). There are 
also other types of BSS introduced in this paper, but they 
rather represent changes of types of BSS connected with the 
dynamical interactions, not new channels of formation. 

Some conclusions in this paper are not general, because 
we discussed here only one test model. In the next papers we 
plan to perform many more simulations with various initial 
parameters (both global and for binaries properties) and we 
will try to study in details the channels of formation of BSS 
and conclusions presented here. 



For EMT channel we noticed groups of BSS which had 
different physical processes responsible for their creation. 
The first group consists of EMT BSS which are created 
through a Roche lobe overflow in compact binaries - short 
period EMT. In this scenario one of the stars in the binary 
starts to leave the main sequence. It increases its radius, 
and at some point, a mass starts to flow to the other star, 
which later on becomes the EMT BSS. The second group 
of EMT BSS is formed in the long period binaries in a dif- 
ferent way. Here, a donor star at some point starts to leave 
the main sequence, enters the giant branch and ejects some 
mass through stellar winds. The star, which later on be- 
comes BSS, gains a mass (about few ~ O.OIM©), which is 
not enough to become a BSS. The donor star then evolves 
to short AGB phase, ejects an envelope. This time, the fu- 
ture BSS, increases the mass by about < O.IMq. The donor 
star quickly becomes a white dwarf, but the star which gains 
mass, in many cases still does not exceeds the turn-off mass. 
This is why these types of EMT are in majority dormant 
BSS. It has to take some time, even several Gyrs, until the 
star exceeds the turn-off mass and becomes a BSS. Thus, 
the most distinct difference, in comparison to short period 
EMT, is that long period EMT gain only small additional 
mass, wait until they exceed the turn-off mass, and in gen- 
eral are short living BSS because they exceed the turn-off 
mass only slightly. EMT are created from around the core 
up to far outside of the half-mass radius. Final EMT po- 
sitions are a little bit more scattered in comparison to the 
initial ones. Some of EMT stayed near and around the core 
and some moved outside the half-mass radius, closer to the 
tidal radius. 

Mass ratios of the long period EMT at the recognition 
time create a distinct trend (see the bottom-right panel in 
Fig. |4]) . The trend is very well reproduced with the Eq. [4] 
The masses of BSS are just slightly larger than the turn-off 
mass, thus in the nominator of Eq. U there is M tU m-of f ■ 
WDs are companions in the long period EMT, so in the 
denominator there ar e calcu lated masses of WDs based on 
IChernoff fc Weinberg! {l990, Tab. 1). The mass ratios for 
long period EMT have a narrow r ange and predictable val - 
ues through the entire simulation. iGeller fc Mathieul (|201ll ) 
reported that BSS in long period binaries in an old (7 Gyr) 
open cluster, NGC 188, have companions with masses of 
about half of the solar mass (surprisingly narrow mass dis- 
tribution). This rules out a collisional origin for these long 
period BSS, because otherwise, for collision hypothesis there 
would be significantly more companions with higher masses. 
It is consistent with a mass transfer origin for the long- 
period blue straggler binaries in NGC 188, in which the 
comp anions would be white d warfs of about half of a solar 
mass l|Geller fc Mathieul 1201 ll ). This finding is very consis- 
tent with our work. We find for long period EMT that their 
mass ratios are in fact narrow. Mass ratios decrease accord- 
ing to the Eq. [3] There are indeed WDs in long period EMT 
with masses of about 0.5M© at the time 7 Gyrs. 

EM are created from compact binaries with eccentrici- 
ties close to zero. For EM BSS we noticed two different for- 
mation processes as well. The first one creates EM through 
the Roche lobe overflow during a semi-detached phase of a 
binary. The second one creates EM BSS from a magnetic 
braking. Although both groups overlaps each other during 
the simulation, EM created from the Roche lobe overflow 
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are created mainly in the beginning of the simulation, and 
EM created from magnetic braking are formed mainly after 
several Gyrs. All EM were created around the core and near 
the half-mass radius. The final positions of EM were more or 
less the same as initial ones with the exception that the final 
positions are a little bit more scattered in the cluster: from 
the core up to even the tidal radius - similarly to EMT. 

Collisional BSS (CBS and CBB) are created in strong 
dynamical interactions. Before the core starts to collapse 
(< 6 Gyrs), CBS and CBB are created in the core from bi- 
naries with mostly large orbital periods (about 10 4 days). 
During and after the core collapse, when the core is get- 
ting denser, the number of CBS and CBB increases. Also, 
they are being created from more compact binaries which 
have orbital periods of about a few days. BSS from dynam- 
ical channels are formed mainly in the core or close to it. 
It is actually not surprising because in order to increase the 
chance of a collision in the dynamical interactions, a dense 
environment is needed. Some of BSS in the dynamical in- 
teractions gained more additional kinetic energy, extended 
their orbits, moved near the tidal radius and beyond. Many 
of them escaped from the star cluster as BSS. 

EMT BSS are the most active in the beginning of the 
star cluster evolution, its population continuously drops but 
nevertheless stays active during the whole cluster evolution. 
The population of EM BSS is not significant in the begin- 
ning. Its peak value in our simulation is about ~ 3 Gyrs, 
when both scenarios work most efficiently for EM - the 
Roche lobe overflow and magnetic braking. Additionally, 
the magnetic braking starts to work for both components 
in binaries (t he turn-off mass dro ps at that time to about 
1.25 M , see iHurlev et all (|2002T l for details). CBS+CBB 
channels starts to play a significant role in the overall pop- 
ulation of BSS i n the post collapse phase. According to 
iFusi Pecci et all l| 19921 ) for less dense GCs, BSS could form 
as evolutionary mergers, in contrary to dense GCs, for which 
BSS could form from dynamical interactions. Our test sim- 
ulation is consistent with it. Before the core collapse, there 
are more evolutionary mergers and the dynamical interac- 
tions become important after the core collapse. Nevertheless, 
we have BSS from different channels during the whole sim- 
ulation. Moreover, we think that the number of BSS from 
evolutionary channels (EM, EMT) depend rather on initial 
binary properties. We expect to have less binaries for EMT 
and EM channels for star clusters where there are less com- 
pact primordial binaries. For simulations with larger initial 
concentrations we also expect to have more BSS from CBS, 
CBB channels and maybe also from CSS. Nevertheless, we 
expect to have BSS from all channels of formation for various 
densities of the star clusters. The only exception could be 
for the CSS channel. In our test model we had only one BSS 
created by the collision of two field stars, and even after the 
core collapse, the number of CSS BSS did not increase. Per- 
haps for star clusters with higher densities CSS will increase 
its overall significance. BSS channels o f formation work si- 
multaneously in di fferent radial parts |Ferraro et alj fl997l : 
iMapelli et all [20061. Indeed, based on our test model, one 
can see that EM and EMT are created more or less in the 
same radial distances from the center. Only BSS from CBS 
and CBB channels are created mainly inside and around the 
core, so in deeper regions than EM and EMT. 

Number of BSS does not correlate with the predicted 



collision rates |Piotto et all |2004|; iLeigh et alj l2008al lbh. 
which is one of the reasons why iKnigge et al. r<|2009h favor 
the mass transfer mechanism as a more important scenario 
for creation of BSS. Our test simulation shows that there 
are indeed many EMT created (49% of all BSS are EMT) 
and that this channel is active for the whole simulation. 
However, the number of EMT drops from the beginning 
and at some point the number of BSS created by merg- 
ers and collisions (EM, CBS, CBB) outnumbers the EMT 
BSS. Moreover, we think that the population of evolution- 
ary BSS depends mainly on initial binary properties and 
thus, for different initial conditions EMT could be indeed 
the most important channel of formation. A smaller num- 
ber of primordial compact binaries could create less EM. A 
stronger dependence on the initial binary properties, rather 
than star cluster global parameters, like concentration, are 
also consisten t with the calculated specific frequencies from 
iFerraro et ah! (120031 ). They found that the largest BSS spe- 
cific frequencies are observed for the star clusters with the 
lowest central density (NGC 288) and the highest central 
density (M80). It suggests, that the number of BSS depends 
rather on the initial binary prop erties, instead of the global 
parameter - the central density. ISollima et alj (|2008h found 
the linear correlation between the number of BSS and the bi- 
nary fraction, which in turn depends mainly on the initial bi- 
nary properties of a star cluster. Thus, it is very important to 
find an observational distinction between the mass transfer 
BS S and the collisiona l ones (fi rst attempts are alrea dy done 
bv lFerraro fc Lanzonil (|2009h V ISollima et al l |2008l ) suggest 
that the strong correlation between the number of BSS and 
the binary fraction is a result of formation channel of BSS as 
the unperturbed evolution of the primordial binaries. They 
did not found correlations with the central density, the con- 
centration, the stellar collision rate, and the half-mass re- 
laxati on time of star clusters. However, later IKnigge et al.1 
(2009) claimed that if one restricts a sample of star clusters 



to these with dense cores, a significant correlation appears 
between the number of BSS and the collision rate. We plan 
to check both, different initial binary properties and differ- 
ent initial global properties, and perform many simulations 
to study how much the population of BSS depends on initial 
binary properties, collision rates and concentrations. 

BSS created as EM or EMT very rarely change their 
types. Only in a few cases EM got into a binary in some 
dynamical interactions and some EMT were disrupted or 
collided with other stars during the dynamical interactions. 
EM do not change types most likely because they are sin- 
gle stars, and the probabilities of the dynamical interactions 
for them is lower than for EMT. Additionally, EM are not 
the most massive objects and the number of EM is small in 
comparison to the number of e.g. WDs at the later stages of 
the star cluster evolution. EMT do not change types, prob- 
ably because of the fact, that most of them are hard enough 
to survive dynamical interactions. Additionally, the number 
of EMT is even smaller than the number of EM, when the 
dynamical interactions become important (> 8 Gyrs). One 
of the subgroups of EMT are long period BSS which should 
favor the higher probabilities of the dynamical interactions. 
However, long period EMT are short living BSS and thus, 
they have small chances to be changed by the dynamical 
interactions. As expected, CBS and CBB channels change 
their types quite often. They are disrupted (if BSS was in 
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a binary) or exchanged with other stars in dynamical inter- 
actions. In many cases, CBS and CBB, stay in binaries or 
get into heavier binaries quickly and thus, they have larger 
masses in general. Additionally, CBS and CBB are created 
mainly in the core or around it. Due to the mass segregation 
they sink to the center of a cluster. All these facts are the 
reasons why the probabilities of the dynamical interactions 
for CBS and CBB are in general higher than for EM or EMT 
and thus they can change types more frequently. 

Many of EM and EMT BSS, which are created after the 
core collapse, had some strong dynamical interactions before 
they became BSS. The strong dynamical interactions means 
that the semi-major axis or the eccentricity was changed by 
some dynamical interaction at least by 10% (arbitrary cho- 
sen value for our test simulation) . It suggests that such EM 
and EMT could be created by possible induced mass trans- 
fers or possible induced mergers. In other words, the cre- 
ation of EM and EMT could be triggered or made possible 
by the dynamical interactions. However, for the simulation 
without the Fewbody, and for simple population synthesis, 
the overall number of EM and EMT were almost the same 
(see Sec. 14.1.71 Fig. Hip . This, in turn, supports the the- 
ory that the dynamical interactions are not important for 
the populations of EM and EMT. On the other hand, for 
a simulation with higher r pmax (higher impact parameters) 
in the MOCCA code, there are significantly more fly-by in- 
teract ions for bina r ies (fo r details see the next paper in the 
series Gicrsz et all (|201ll . Fig. 7)). In such a case there were 
also created many more BSS. This, in turn, suggests that the 
dynamical interactions could have some influence on the cre- 
ation of BSS. We plan to check in the next papers whether 
the dynamical interactions have indeed some implications 
for the population of EM or EMT. Maybe the star clusters 
with higher concentrations, or for which dynamical inter- 
actions are more important in overall evolution, will have 
larger population of EM and EMT. 

A bimodal distribution of BSS is present i n ma n y star 
clusters. It was discovered b y Ferraro et alj (1 19931 . 19971 ) 
for M3 and by IZaggia et all l| 19971 ) for M55. Later, the 
bimodality was observed for many other clusters and for 
a few of them it was not present, e.g. for NGC 2419 
l|Dalessandro et al.ll200sl ; IContreras Ramos et al.ll2012l ). We 
were also able to see weak signs of the bimodal distribution 
of BSS in our test si mulation. The bimod al distribution of 
BSS was explained bv lFerraro et al.l (|l997l ) as a result of dif- 
ferent processes forming BSS in different parts of the cluster 
- mass transfers for the ou ter BSS and s t ellar collision s for 
BSS in the center. Later, iMapelli et ail (|2004 l200rj ) and 
lLanzoni et all ll2007l) . showed that the bimodal distribution 
in the cluster cannot be explained only by the collisional 
scenario, when BSS are created in the center of the cluster 
and some of them are ejected to the outer parts of the sys- 
tem. This process is believed to be not efficient enough and 
~ 20 — 40% of BSS have to be created in the peripherals in 
order to get the required number of BSS in the star clus- 
ters. It is believed that in the outer parts of the star cluster 
binaries can start the mass transfer in isolation without suf- 
fering from the energetic dynamical interactions with field 
stars. The last explanation is very consistent with our test 
simulation. BSS from EM and EMT channels are created 
from the core up to the outside of the half-mass radius, and 
CBS+CBB are created mainly inside and around the core 



(see Sec. 14.2. . Some of CBS and CBB are indeed ejected to 
more extended orbits. In total, the number of BSS outside 
the half-mass radius consist of all channels (EM, EMT and 
CBS+CBB). Moreover, the bimodality of BSS is not present 
always in our test simulation, but it seems to form after the 
core collapse and becomes visible after that. If this could 
be confirmed in our next planned simulations with differ- 
ent initial conditions, the bimodality could be a tracer of a 
dynamical status of a cluster. 

We investigated lifetimes of BSS and noticed that for 
some BSS there is a significant delay between the creation 
time (the time of the last merger or the last mass trans- 
fer) and the recognition time (the time when a star actually 
exceeded the turn-off mass). For some BSS this delay can 
last even several Gyrs. Such a delay was unexpected to find. 
BSS, for which such a delay exists, we call dormant BSS. 
The number of dormant EM and EMT is significant. For all 
149 EM BSS there were 45 dormant EM BSS (30%), for all 
231 EMT BSS there were 60 dormant EMT (26%) and for 
all 95 CBS and CBB only 7 were dormant (7%). For the 
total 476 BSS there were overall 112 dormant BSS which is 
24%. 

There is a number of BSS which escaped from the 
star cluster. For EM and EMT, this process seems to be 
not important because only respectively 3% and 9% es- 
caped as BSS. However, for CBS and CBB channels for the 
whole simulation, 40 BSS (from the total 95) escaped as 
BSS, which gives 43% efficiency. The number of CBS and 
CBB escapers increase even more if one narrow the time 
to the post collapse. After the core collapse, 60% of CBS 
and CBB escaped as BSS from the star cluster. There are 
known BSS w hich are found in the halo and in the bulge 
of the Galaxy jlBragaglia et al.ll2005l ; iFuhrmann et al.ll201ll ; 



IClarkson et al.l 2011 ), an d there are also known fast moving 
BSS (|Tillich et al.ll201ol ). These BSS were probably created 
from CBS or CBB channels, because from the dynamical in- 
teractions these stars could get high escape velocities. Our 
test model seems to show that CBS and CBB escapers could 
be important. Also, it seems that for BSS escapers channel 
of formation determines what is the cause that BSS leaves 
the system. If EM and EMT are ejected from the star clus- 
ter, it is because of the relaxation processes. CBS and CBB, 
in the contrary, leave the system because of the dynami- 
cal interactions - they gain additional kinetic energy. CBS 
and CSS are fast escapers, whereas EM and EMT are slow 
ones. For EM and EMT there is significant delay between 
the recognition time and the escape time (see Fig. I16p . We 
plan to check in the next paper how significant BSS esca- 
pers will be for all channels of formation for different initial 
parameters. 

5.1 Future plans 

In the next papers we plan to perform a systematic study for 
different binary properties, different IMF, concentrations, 
metallicities and many other star cluster parameters. We 
plan to check how it can influence the overall population of 
BSS or the population of the different channels of formation. 
We plan to search for correlations between results of the 
simulations and observational data. 

The Fewbody code allows to deal with any kind of hier- 
archies like triples, quadruples and higher hierarchies. How- 
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ever, currently the MOCCA code is unable to deal with such 
type of objects. It works only with single and binary stars. 
One of the significant new features would be to implement 
procedures to deal with the higher hierarchies. 

Another major new feature would be parallelization of 
the source code. The usage of the Fewbody slowed down 
the code. In order to regain, as much as possible, the pre- 
vious performance of the MOCCA code, the next step is 
to parallelize the main bottlenecks of the code. Executing 
the dynamical interactions in parallel could give some speed 
up. A natural choice is to use OpenMP because the main 
procedures are just loops which sequentially execute specific 
functions. 
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